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1  Motivation  and  Objectives 

One  of  the  most  promising  techniques  for  the  prediction  of  turbulent  flows  is  that  of  Large  Eddy 
Simulation  (LES),  in  which  an  under-resolved  representation  of  the  turbulence  is  simulated  numer¬ 
ically  by  modeling  the  effects  of  the  unresolved  small-scales  on  the  simulation.  Such  simulations 
have  been  applied  in  several  flows  with  reasonable  success.  However,  there  are  several  outstand¬ 
ing  problems  that  need  to  be  addressed  before  LES  can  fulfill  its  promise  as  a  tool  for  turbulence 
prediction  in  engineering  flows.  The  most  serious  problems  limiting  the  usefulness  of  LES  is  the 
representation  of  turbulence  near  walls  and  other  strong  inhomogeneities  and  the  dependence  of 
models  on  the  filter  and/or  numerical  discretization.  Also  important  are  the  treatment  of  inhomo¬ 
geneous  filters  and  the  lack  of  understanding  of  the  modeling  errors  and  their  impact. 

The  optimal  LES  formulation  [9,  10]  provides  a  rigorous  framework  in  which  to  address  these 
issues  and  to  develop  and  analyze  LES  models  and  simulations.  The  objective  of  this  research  is  to 
develop  new  LES  modeling  techniques  to  address  the  difficulties  cited  above  by  using  the  optimal 
LES  formulation. 


2  Background  &  Approach 


The  starting  point  for  the  development  of  LES  is  the  definition  of  a  spatial  filter  r,  which  can  be 
applied  to  the  Navier-stokes  equations  to  obtain  an  equation  for  the  filtered  velocity  uc. 

duj  =  dujUj  dp  1  d2Uj 

dt  dxj  dxi  R  edxjdxj 

Where  Mj  is  the  sub-grid  model  (force)  term,  which  includes  the  divergence  of  the  sub-grid  stress 
as  well  as  terms  that  arise  when  the  filter  does  not  commute  with  differentiation.  The  problem  in 
LES  of  course  is  to  model  Mj.  A  very  important  result  of  our  research  [9]  is  that  an  LES  w  will 
match  the  one-time  statistics  of  filtered  turbulence  u  if  and  only  if  the  model  rrii(w )  of  Mj  is  given 
by 

rrii(w)  =  {Mi(u)\u  =  w)  (2) 
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This  model  also  minimizes  the  difference  between  Mi  and  m,  (in  the  mean-square  sense),  and  so 
this  model  has  all  the  properties  that  one  could  ask  of  a  sub-grid  model.  We  therefore  call  it  the 
ideal  sub-grid  model. 

Unfortunately,  the  conditional  average  in  (2)  cannot  practically  be  determined,  since  the  conditions 
are  that  the  entire  filtered  velocity  field  match  the  entire  LES  field.  However,  it  can  be  estimated 
using  stochastic  estimation  [2]  which  is  a  well-established  technique  for  estimating  conditional 
averages.  The  result  is  a  class  of  estimation  based  LES  models  as  first  proposed  by  Adrian  [1] 
To  perform  stochastic  estimation  one  requires  large  quantities  of  two-point  correlation  data.  For 
example,  in  linear  estimation  one  needs  (uj{pd)uk{x!'))  ar*d  (Mj(x)u*(x")). 

To  obtain  the  required  data,  we  resort  to  direct  numerical  simulation  (DNS),  which  is  limited  to  low 
Reynolds  numbers,  and  experimental  velocity  field  measurements,  which  are  currently  limited  to 
two-dimensional  measurements.  By  using  reliable  empirical  data  for  the  required  correlations,  the 
optimal  modeling  approach  can  be  developed,  tested  and  validated  in  the  absence  of  uncertainties 
associated  with  modeling  the  correlations  theoretically.  Once  the  formulations  has  been  developed 
and  validated,  truly  predictive,  generally  applicable  models  can  be  developed  based  on  theoretical 
approximations  to  the  required  correlations. 


3  Supported  Research 


To  pursue  the  objectives  defined  above,  a  number  of  research  activities  were  pursued  under  the 
current  grant.  These  are  described  briefly  below  and  in  more  detail  in  the  following  subsections, 
and  the  referenced  publications. 


1.  Correlation  Data  Acquisition:  To  develop  and  test  the  optimal  LES  formulation,  exten¬ 
sive  multi-point  correlation  information  was  needed.  Since  the  wall-modeling  problem  is 
of  prime  importance,  we  concentrate  on  near-wall  correlations  from  turbulent  channel  flow. 
Data  was  generated  both  experimentally,  using  a  novel  PIV  technique  [6],  and  numerically, 
using  direct  numerical  simulation  (DNS).[8]  This  data  is  also  being  used  to  formulate  and 
check  theoretical  descriptions  of  the  two-point  correlation  (see  item  5  below). 

2.  Development  and  Testing  of  Optimal  LES  Formulations:  It  had  previously  been  shown 
that  optimal  LES  based  on  Fourier  cut-off  filters  performed  very  well.  But,  Fourier  cutoff 
filters  are  not  practical  for  use  in  complex  geometries.  So,  an  optimal  LES  formulation 
based  on  finite  volume  filters  was  developed  and  tested  (see  [16]  and  Appendix  A).  Further, 
for  several  reasons,  it  would  be  advantageous  to  develop  optimal  models  based  on  finite 
difference  discretizations,  rather  than  finite  volume,  and  such  formulations  have  also  been 
developed  and  tested. 

3.  Filtered  Wall  Treatment:  One  of  the  difficulties  associated  with  wall-bounded  flows  is  that 
if  one  sharply  defines  the  wall  and  the  no-slip  boundary  condition,  it  requires  that  the  filter  be 
extremely  inhomogeneous  near  the  wall.  This  causes  a  number  of  problems.  Alternatively, 
the  boundary  can  itself  be  filtered,  producing  the  LES  analog  of  an  embedded  boundary 
method.  Wall  models  for  this  approach  have  been  developed  using  optimization  techniques. 
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4.  Theoretical  Optimal  LES  for  Isotropic  Turbulence:  To  eliminate  the  need  for  extensive 
empirical  input  to  optimal  modeling,  the  required  correlations  need  to  be  determined  theoret¬ 
ically.  Using  a  combination  of  Kolmogorov  scaling,  small-scale  isotropy,  the  quasi-normal 
approximation  and  a  dynamic  approach,  theoretical  optimal  LES  models  for  the  finite  vol¬ 
ume  formulation  were  developed.  The  resulting  models  are  valid  provided  the  Reynolds 
number  is  sufficiently  large,  and  small-scale  isotropy  is  a  valid  approximation.  These  mod¬ 
els  have  also  been  tested. 

5.  Theoretical  Correlations  for  Wall-Bounded  Turbulence:  Near  walls,  the  assumption  of 
small-scale  isotropy  and  filter  widths  in  the  inertial  range  are  invalid.  A  representation  for 
the  inhomogeneity  and  anisotropy  of  the  multi-point  correlations  in  wall-bounded  turbulence 
(particularly  in  the  log-layer)  is  needed  for  use  in  the  development  of  optimal  models.  Such 
a  representation  is  being  pursued  based  on  the  developments  by  Procaccia’s  group[4, 3, 11], 
the  similarity  forms  of  Oberlack,[13]  consistency  with  the  Navier-Stokes  equations  and  the 
quasi-normal  approximation.  Further,  the  DNS  and  experimental  data  described  in  (1)  are 
being  used  for  guidance. 


3.1  Correlation  Data  Acquisition 


To  acquire  the  correlation  data  needed  to  address  the  modeling  of  wall-bounded  turbulence,  a 
new  DNS  of  turbulent  channel  flow,  and  a  laboratory  experiment  on  turbulent  channel  flow  were 
conducted. 


3.1.1  Channel  Flow  DNS 


A  turbulent  channel  flow  at  Rer  =  940  in  a  very  large  spatial  domain  (Lx  =  Snh,  Lz  =  3.257 xh, 
where  h  is  the  channel  half-width)  was  conducted.  The  large  spatial  domain  allows  more  reliable 
statistics  to  be  gathered,  and  facilitates  the  study  of  large  length-scale  phenomena  known  to  be 
present  in  the  near-wall  region.  Further,  this  large-domain  channel  provides  a  valuable  validation 
target  for  wall-bounded  turbulence  LES.  In  particular  the  capture  of  long  wavelength  phenomena 
by  the  LES  can  be  assessed. 

The  simulations  were  performed  at  the  San  Diego  Supercomputer  facility  and  at  the  DOD  HPC  fa¬ 
cility  at  NAVO.  The  simulation  produced  a  library  of  approximately  100  DNS  fields  which  are  cur¬ 
rently  being  post-processed  for  a  number  of  purposes,  including  the  correlation  representation  dis¬ 
cussed  in  section  3.5.  Furthermore  the  data  is  in  the  process  of  being  made  available  to  researchers 
in  the  turbulence  community  through  the  website  http:  /  /da vinci .  tam.uiuc . edu.  Fi¬ 
nally,  the  first  paper  published  on  the  data  has  appeared  [8]. 


3.1.2  Channel  Flow  Experiments 


To  support  the  development  of  optimal  LES  models,  detailed  PIV  measurements  of  turbulent  chan¬ 
nel  flow  were  made.  In  addition  to  standard  velocity  measurements,  which  are  used  to  determine 
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the  multi-point  velocity  correlations  required  in  the  formulation  velocity  time-derivative  measure¬ 
ments  were  made  using  a  novel  PIV  technique,  and  correlations  of  the  time-derivative  were  com¬ 
puted  also  to  support  modeling  of  LES  dynamics.  The  experimental  technique  is  described  briefly 
below,  and  in  more  detail  in  [6], 

Measurements  of  instantaneous  time-  and  bulk-convective-derivative  fields  have  been  made  over 
a  broad  range  of  Reynolds  numbers  to  support  the  development  of  optimal  LES  models.  A  two- 
CCD-camera  PIV  arrangement  was  designed  to  acquire  large  ensembles  of  time-resolved  velocity 
data  at  moderate-to-high  Reynolds  numbers  in  the  streamwise-wall-normal  plane  of  fully  devel¬ 
oped  turbulent  channel  flow.  The  cameras  are  focused  on  the  same  field  of  view  (temporal  deriva¬ 
tive)  or  shifted  in  the  streamwise  direction  (bulk  convective  derivative),  and  the  acquisition  of  the 
second  is  delayed  in  time  by  r  (typically  a  fraction  of  the  Kolmogorov  time  scale,  4),  to  allow  the 
computation  of  the  associated  time  derivative.  Polarization  filtering  is  used  to  separate  the  particle 
images  viewed  by  each  camera.  Rigorous  techniques  have  been  developed  to  compensate  for  the 
noise  inherent  in  the  PIV  velocity  measurements. 

Recent  effort  has  focused  on  the  bulk  convective  derivative  because  comparison  of  the  streamwise 
spectra  of  dui/dt  and  Dui/Dt  indicates  that,  at  all  scales,  the  temporal  derivative  of  velocity 
is  dominated  by  bulk  convection,  not  evolution.  These  convection  effects  mask  the  underlying 
dynamics  of  interest  in  optimal  LES  development. 

Reynolds-number  scaling  of  the  bulk  convective  derivative  statistics  has  been  assessed  for  both  un¬ 
filtered  and  filtered  data.  The  intent  of  this  analysis  is  to  determine  the  most  appropriate  time  scale 
for  achieving  Reynolds-number  similarity  in  the  bulk  convective  derivative  statistics.  Figure  1(a) 
shows  the  root-mean-square  of  Dbv/Dt  (unfiltered)  scaled  with  the  friction  velocity  («*)  and  the 
outer  time  scale  ( t„  =  u*/h)  as  a  function  of  wall-normal  position  (Dbv/Dt  has  the  strongest 
Reynolds  number  dependence).  Near  the  wall,  this  scaling  does  not  remove  the  strong  Reynolds- 
number  dependence  in  the  RMS,  but  further  from  the  wall  ( y/h  >  0.4),  this  scaling  works  well. 
This  is  the  best  scaling  found  for  this  quantity.  In  figure  1(b)  the  RMS  time  derivative  of  the  filtered 
fields  is  shown  (top-hat  filter  of  width  A  =  0.2 h).  The  strong  Reynolds-number  dependence  noted 
in  the  unfiltered  profiles  is  now  absent.  Scaling  of  the  two-point  correlation  functions  of  DbUi/Dt 
yields  similar  results:  unfiltered  correlations  show  Reynolds-number  dependence  when  scaled  with 
t0,  while  filtering  removes  this  dependence,  yielding  a  consistent  scaling  of  the  correlations.  It  is 
the  time  derivative  of  the  filtered  velocity  that  is  of  interest  in  LES,  so  the  scalings  based  upon 
tQ  are  appropriate  for  the  optimal  LES  formulation.  This  information  will  be  used  to  develop  a 
comprehensive  Reynolds-number-scaling  framework  for  extending  optimal  formulations  to  higher 
Reynolds  numbers. 


3.2  Development  and  Testing  of  Optimal  LES  Formulations 

A  finite  volume  formulation  of  optimal  LES  was  developed  previously  [10],  but  had  not  been 
tested.  This  formulations  has  been  tested  and  refined.  Also,  it  was  recognized  that  a  finite  differ¬ 
ence  formulation  of  optimal  LES  would  be  advantageous,  and  this  has  also  been  developed  and 
tested,  though  more  work  on  refinement  remains  to  be  done.  These  activities  are  discussed  below. 
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Figure  1:  Scaling  of  T)bv/DtRMs  profiles  with  u*  and  ta.  (a)  Unfiltered,  (b)  Filtered. 


Figure  2:  Three-dimensional  energy  spectrum  E(k),  filtered  DNS  compared  with  optimal  LES  with 
2-cell  (S2)  and  4-cell  (S4)  stencils  and  dynamic  Smagorinsky. 


3.2.1  Finite  Volume  Optimal  LES 


In  the  finite  volume  LES  formulation  of  Langford  [10],  the  “filter”  is  defined  to  be  the  average 
over  discrete  volumes.  As  in  a  standard  finite  volume  numerical  discretization,  the  evolution  of  the 
volume  averaged  velocity  requires  the  evaluation  of  momentum  fluxes  through  the  boundaries  of 
each  discrete  volume.  In  the  LES  context,  these  cannot  be  estimate  using  standard  finite  volume 
reconstruction  techniques  because  the  cell  size  (filter  width)  is  not  small  compared  to  the  smallest 
scales  over  which  the  velocity  varies  (the  Kolmogorov  scale).  In  short,  the  velocity  is  not  smooth 
on  the  scale  of  the  volumes.  Instead,  we  use  the  optimal  LES  technique  to  estimate  the  fluxes. 
Thus,  the  usual  convergence  considerations  normally  used  to  develop  numerical  methods  are  re¬ 
placed  with  the  requirement  that  the  flux  estimates  be  consistent  with  the  statistical  properties  of 
turbulence. 

This  finite  volume  LES  formulation  was  tested  in  isotropic  turbulence  using  DNS  data  to  determine 
the  required  correlations.  The  results  are  quite  good  as  shown  by  the  spectra  plotted  in  figure  3.2.1. 
It  was  found  that  a  simple  stencil  in  which  two  cells  on  either  side  of  a  face  are  used  to  estimate 
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the  flux  is  adequate  to  yield  very  good  results.  It  was  also  found  that  a  staggered  grid  formulation 
produces  a  much  better  LES.  The  details  of  this  study  are  provided  in  [16],  a  preprint  of  which  is 
included  as  Appendix  A. 


3.2.2  Finite  Difference  Optimal  LES 

One  issue  which  must  always  be  addressed  in  LES  is  that  it  is  generally  not  the  filtered  turbulence 
that  is  of  interest,  but  rather  the  real  unfiltered  turbulence.  Thus,  when  interpreting  LES  results,  it 
is  necessary  to  account  for  the  missing  small-scales.  For  example,  the  turbulent  kinetic  energy  will 
be  missing  the  contribution  from  the  unresolved  small  scales. 

If  one  is  interested  in  statistical  correlations  of  the  velocity,  then  an  LES  definition  that  preserves 
the  statistics  of  the  velocity  would  be  particularly  useful,  since  it  would  avoid  these  difficulties. 
One  such  LES  filter  definition  is  just  a  sampling  of  the  turbulence  on  a  grid  of  points.  This  is  a 
non-invertible  linear  transformation  that  can  be  used  as  an  LES  “filter,”  though  it  does  not  have 
the  usual  filter  properties  of  eliminating  small-scale  fluctuations.  Indeed,  such  a  mapping  can  be 
considered  to  exhibit  aliasing  of  the  small  scales.  None-the-less,  the  ability  to  preserve  velocity 
statistics  in  this  formulation  make  it  an  attractive  target  for  LES  development. 

The  evolution  equations  for  the  point-sampled  velocity  are  given  by 


dv?i 

duiuk 

j 

dp 

j 

+  u 

d2Ui 

dt 

dxk 

dxi 

dxkdxk 

where  j  is  the  index  for  the  grid  of  points.  Note  that  the  terms  on  the  right  hand  side  all  involve 
spatial  derivatives,  which  cannot  be  directly  evaluated  from  the  point  sampled  fields.  Because 
the  grid  spacing  is  assumed  to  be  large  compared  to  the  smallest  turbulence  scale,  the  usual  finite 
difference  approximations  of  the  derivative  are  invalid.  Instead,  we  use  optimal  LES  modeling  to 
estimate  these  derivatives.  This  is  very  like  the  finite  volume  formulation,  though  we  are  estimating 
derivatives  in  this  case,  rather  than  fluxes. 

There  were  several  subtleties  that  needed  to  be  addressed  in  developing  the  optimal  estimates  of 
the  derivatives.  First  among  these  is  that  the  nonlinear  term  is  energy  conserving  in  this  formu¬ 
lation.  All  energy  dissipation  occurs  through  the  viscous  term.  It  was  necessary  for  stability  and 
robustness  that  the  estimated  nonlinear  term  be  constrained  to  be  energy  conserving  as  well.  This 
was  easily  accommodated  as  part  of  the  estimation  procedure.  Another  issue  is  the  enforcement 
of  continuity.  The  presence  of  significant  aliasing  in  the  LES  filter  definition  makes  continuity  a 
very  weak  constraint  on  the  point  sampled  velocity.  Thus,  imposing  a  continuity  constraint  on  the 
sampled  velocity  is  apparently  not  valid.  But  the  pressure  must  still  be  determined.  For  our  prelim¬ 
inary  studies,  an  approximate  continuity  constraint  was  applied  and  the  pressure  was  determined 
from  this  constraint.  Better  approaches  are  being  explored. 

A  preliminary  result  for  the  spectrum  of  the  sampled  velocity  is  shown  in  figure  3,  along  with 
that  for  the  sampled  DNS  of  isotropic  turbulence.  The  results  are  not  nearly  as  impressive  as  the 
finite  volume  results  discussed  above,  but  at  least  part  of  this  can  be  attributed  to  the  treatment  of 
continuity.  The  finite  difference  optimal  LES  approach  is  continuing  to  be  refined. 
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Figure  3:  Three  dimensional  energy  spectrum  in  forced  isotropic  turbulence  with  a  point-sampled 
LES  filter.  Shown  are  filtered  DNS  and  LES  results. 

3.3  Filtered  Wall  Treatment 

It  was  pointed  out  by  Das  &  Moser[7]  that  the  wall  of  a  wall  bounded  flow  could  be  filtered  to 
avoid  the  use  of  strongly  inhomogeneous  filters  near  the  wall.  In  this  approach,  a  homogeneous 
(or  nearly  homogeneous)  filter  is  applied  to  an  extended  domain  including  a  region  exterior  to  the 
wall  in  which  u  =  0.  When  such  a  filter  is  applied  to  the  Navier-Stokes  equations  (extended  in  to 
the  interior),  the  result  is 

dui  -  dui_  _  dp  1  d2Uj  dTjj 

dt  dxj  dxi  Redxjdxj  *  dxj 

where  bi  is  the  boundary  term  arising  from  the  filtering  of  the  boundary  and  is  the  usual  subgrid 
stress  term.  The  boundary  term  is  expressed  as 

&i(x)  =  f  <Tij(x')njG(x  —  x')  dx1 
j  dR 

where  a  is  the  stress  at  the  boundary,  including  pressure  and  viscous  stress,  dR  is  the  boundary 
of  the  fluid  region  R  and  rij  is  the  unit  normal  to  the  boundary.  To  determine  the  stresses,  we 
note  that  the  velocities  should  remain  zero  exterior  to  the  flow  domain.  The  stresses  can  thus  be 
determined  to  minimize  the  velocities  in  the  exterior.  In  particular,  we  minimize  the  cost  function: 

E  —  J  \v,i\2  +  a  dy,  where  a  is  a  numerical  constant  selected  to  be  of  order  At2. 

To  test  this  approach,  a  filtered  boundary  simulation  of  the  turbulence  in  a  channel  at  ReT  —  180 
was  performed,  with  no  filtering  in  the  horizontal  directions,  and  a  Fourier  cutoff  filter  with  kc  = 
371/ <5  in  the  wall-normal  direction  (A y+  =  1.5).  This  is  a  sufficiently  fine  filter  that  the  primary 
effect  is  to  smear  the  boundaries,  so  it  could  be  considered  to  be  a  filtered-boundary  DNS.  The 
results  of  this  simulation  are  compared  to  those  of  [12]  in  figure  4,  and  the  agreement  is  excellent. 

These  results  indicate  that  the  filtering  and  wall-modeling  approach  are  viable,  at  least  in  the  con¬ 
text  of  a  well  resolved  simulation.  However,  the  real  application  for  this  approach  is  in  LES,  where 
the  filter  width  will  be  large  in  wall  units.  In  this  case,  the  wall  model  must  represent  signifi¬ 
cantly  more  of  the  near-wall  physics,  particularly  the  near-wall  production  process.  To  evaluate 
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Figure  4:  Mean  velocity  (left)  and  rms  velocities  (right)  in  a  turbulent  channel  at  ReT  =  180; 

- ,  DNS  of  [12]; - filtered  boundary  DNS  using  a  Fourier  cutoff  filter  with  kc  =  371  /S 

(A y+  =  1.5).  On  the  right,  streamwise  (top  curves),  spanwise  (middle  curves)  and  wall-normal 
(lower  curves)  rms  velocities  are  shown. 


Mean  velocity  RMS  uf 


Figure  5:  Mean  velocity  and  streamwise  rms  velocity  for  LES  of  channel  flow  at  Rer  =  590.  with 
Fourier  cutoff  filters  in  a  spatial  directions,  thus  the  walls  are  also  filtered.  Shown  are  the  DNS 
data,  the  filtered  DNS  and  the  LES.  Filter  widths  were:  A+  =  116,  Ay  =  37+  and  A+  =  58. 
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the  ability  of  the  model  described  above  to  represent  the  near-wall  process,  an  channel  flow  LES 
was  formulated  using  a  Fourier  cutoff  filter  in  the  wall  normal  direction,  in  addition  to  the  ho¬ 
mogeneous  directions.  In  plus  units  the  filter  widths  in  the  streamwise  wall-normal  and  spanwise 
directions  were  A+  =  116,  Ay  =  37+  and  A+  =  58.  The  Reynolds  number  was  ReT  =  590, 
the  same  case  investigated  by  Voelker  et  a/[15].  An  optimal  LES  model  based  on  the  DNS  data 
of  [12]  was  used  for  the  interior,  with  the  same  formulation  suggested  by  Voelker  et  al.  The  wall 
treatment  discussed  above  was  also  used.  The  results  are  shown  in  figure  5.  The  mean  velocity  is  in 
only  fair  agreement  with  the  filtered  DNS  data,  and  the  rms  streamwise  velocity  is  systematically 
over-predicted,  particularly  near  the  wall. 

It  appears  that  these  shortcomings  arise  because  the  wall  model  is  optimized  to  minimize  “leakage” 
of  energy  and  momentum  into  the  exterior,  rather  than  the  transport  of  energy  (or  Reynolds  stress) 
into  the  fluid  domain,  as  suggested  by  the  results  of  Voelker.  Both  these  requirements  are  important, 
so  a  hybrid  approach  in  which  both  are  minimized  simultaneously  is  being  developed. 


3.4  Theoretical  Optimal  LES  for  Isotropic  Turbulence 

We  seek  to  develop  the  correlations  required  to  do  finite  volume  LES  analytically  based  on  the 
assumption  of  a  Kolmogorov  inertial  range  and  standard  two-point  modeling  assumptions,  partic¬ 
ularly  the  quasi-normal  approximation.  The  finite  volume  optimal  LES  is  selected  as  a  target  here 
because  the  goal  is  a  generally  applicable  modeling  approach,  which  makes  spectral  representa¬ 
tions  and  filters  inappropriate. 

In  the  finite  volume  optimal  LES  formalism,  we  need  to  estimate  the  integrated  fluxes  through  a 
surface  in  terms  of  the  velocities  averaged  over  discrete  volumes,  and  these  estimates  must  be  at 
least  quadratic  to  represent  the  Navier-Stokes-like  terms  of  the  LES.  Thus,  our  estimate  is  of  the 
form: 

I  «i(x)u,(x)  dx  =  ^2  Lij(s,  v)  /  ttj(x)  dx  (5) 

Js  V  '*v 

+  H  Qijk(s,v i,t/2)  [  Ujix^dx1  [  uk{x2)dx2  (6) 

vl,v2  Jv'  Jv 2 

where  s  is  a  surface  (volume  face)  and  v  a  volume,  and  us  is  the  (outward  directed)  component  of 
u  normal  to  the  surface  s.  To  determine  the  estimation  coefficients  L  and  Q,  we  solve  a  system  of 
equations  in  terms  of  the  integrated  multi-point  correlations: 

Iu(v',s)  =  ^2Lij(s,v)Iij(v\v)+  23  Qijk{s,Vi,V2)Iijk(v',Vi,V2)  (7) 

V 

V2i  S)  =  Ljj (s,  v)IlmiiVu  v2iV)  ^  Qijk(si  VU  V2)Ilmjk(Vli  V2i  ^2)  (8) 

V  VI  ,V2 

where  the  integrated  correlations  / 1  to  75  are  given  by 


Ili(v',s)  =  J '  £(ui(x')ui(x)us(x))dxdx' 
IijW,v )  =  j(ui{x!)uj{x))dxdx! 
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(11) 


Ifjk(v',v1,v2)  =  f  [  f  (ui(x')uj(x.1)uk(x2))  dx2  dx1  dx' 

Jvf  Jv  1  Jv 2 

lLij(v[,v'2ls)  =  /,/,  ^(^(xl,)wm(x2,)ui(x)ws(x))dxdx2'dx1/  (12) 

^1,^2)  =  Jv>  Jv>  £  £  (u;  (x1,)'Um  (x2/)u,  (x1  )«fc(x2))  dx2  dx1  dx2'  dx1'  (13) 

Thus,  there  are  three  correlation  tensors  needed  to  determine  the  estimation  coefficients: 

-M1*1)  =  (wi(x)wj(x1))  (14) 

Tiifc(r\r2)  =  (ui(x)ui(x1)wfc(x2))  (15) 

■^iw(rl>r2,r3)  =  (wi(x)wi(x1)wfc(x2)^(x3)),  (16) 


where  we  have  used  homogeneity  to  express  the  dependence  of  the  correlation  on  location  differ¬ 
ences,  rl  =  x  —  x*. 


3.4.1  Correlations 


To  determine  these  correlations,  we  can  make  use  of  the  isotropy  of  the  tensors,  and  we  will  assume 
that  all  the  spatial  separations  are  small  enough  to  be  in  the  Kolmogorov  inertial  range  of  an  infinite 
Reynolds  number  isotropic  turbulence.  From  this  assumption,  we  can  determine  the  second  and 
third-order  longitudinal  structure  functions: 

S2(r')  =  ((t1|l(x)-u,(X1))2)  =  Cl£2/3(r1)V3  (17) 

S3(r‘)  =  ((ti||(x)-ti|(x1))3)  =  -|£r1  (18) 

where  r*  =  |r*|  is  the  magnitude  of  the  separation  vector,  u\\  is  the  velocity  component  in  the 
direction  of  the  separation  vector  and  C\  is  the  Kolmogorov  constant. 

The  structure  function  relations  along  with  isotropy  allow  us  to  determine  two  useful  correlations. 
First,  the  most  general  isotropic  two-point  correlation  Rij  that  is  consistent  with  (17)  is  given  by 

Rij( r1)  =  u2Sij  +  9±<?l\rl)-A/\r\r)  -  4(r1)25ii)  (19) 

Similarly,  53  is  related  to  the  third-order  two-point  correlation: 

hjkir1)  =  (wi(x)wi(x)ufe(x1))  =  Tijk(0,  r1).  (20) 

and  the  most  general  isotropic  form  consistent  with  (18)  and  continuity  is: 

Mr1)  =  ^  (W  -  |(far‘  +  V.1))  <21> 

This  is  precisely  the  correlation  needed  for  the  integral  1 1  (9).  However,  the  more  general  third- 
order  three-point  correlation  is  needed  for  J3  (11).  And  the  fourth  order  four-point  correlation  is 
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Figure  6:  Three-dimensional  energy  spectrum  E(k),  filtered  DNS  compared  with  Optimal  LES 
with  4-cell  stencils,  developed  from  theory  and  from  DNS  statistical  data. 


needed  for  J4  (12)  and  75  (13).  Further  modeling  assumptions  will  be  required.  Particularly,  the 
quasi-normal  approximation  can  be  used  to  determine  Fijki  in  terms  of  Rij.  The  result  is: 

Fijki( r\  r2,  r3)  =  ^(r1)#^3  -  r2)  +  #ife(r2)ify(r3  -  r1)  +  Rii(r3)Rjk(r2  -  r1).  (22) 

While  it  is  commonly  used,  the  quasi-normal  approximation  has  had  mixed  success  in  turbulence 
modeling.  Particularly,  when  used  as  a  model  to  close  the  third-order  spectrum  equations,  the 
quasi-normal  approximation  leads  to  an  unrealizable  energy  spectrum.  The  reason  is  that  the  small 
errors  in  the  quasi-normal  approximation  accumulate,  and  realizability  of  the  energy  spectrum  is 
sensitive.  Fortunately,  in  the  current  case,  there  is  no  analog  to  the  realizability  constraint.  Our 
model  is  for  the  filtered  fluctuating  velocity  equation,  which  can  be  used  to  determine  the  spectrum 
directly  in  a  way  that  is  inherently  realizable  (because  it  is  computed  from  the  velocities).  Whether 
the  errors  in  the  quasi-normal  approximation  are  significant  for  the  current  purposes  can  only  be 
determined  by  testing  the  resulting  models. 

The  third-order  three-point  correlation  is  problematic.  It  cannot  be  determined  theoretically  from 
the  considerations  discussed  above.  Instead  we  observe  that  73  (as  well  as  7 2  and  75)  can  be  recast 
as  a  correlation  among  the  LES  state  variable  (volume  averaged  velocities).  Thus  it  is  possible  to 
determine  72,  73  and  75  from  a  running  LES.  This  allows  for  the  possibility  that  these  quantities 
can  be  determined  dynamically  as  a  simulation  is  running.  This  is  the  approach  that  is  currently 
being  pursued  for  73. 

While  this  approach  is  being  refined,  a  test  of  the  modeling  approximations  described  above  was 
conducted  by  computing  73  from  DNS  while  using  the  models  described  above  to  determine  the 
remaining  terms.  The  resulting  model  was  tested  in  isotropic  turbulence  using  the  same  case  and 
same  filter  size  used  in  [16]  (see  Appendix  A).  The  results  are  shown  in  figure  6,  with  comparison 
to  filtered  DNS  and  to  the  DNS-based  model  discussed  in  [16].  The  results  with  this  model  are  not 
quite  as  good  as  the  DNS-based  model,  but  there  is  a  small  inconsistency  associated  with  using  the 
DNS  data  for  73,  and  the  high-Reynolds  number  theory  for  the  rest.  This  may  well  be  responsible 
for  the  this  somewhat  degraded  performance.  None-the-less,  the  results  suggest  that  the  approach 
to  modeling  the  correlations  discussed  here  is  viable. 
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Figure  7:  Contours  of  normalized  error,  4>n,n>  as  a  function  of  the  (left)  streamwise  separations 
and  (right)  spanwise  separations  (along  horizontal  coordinate)  and  wall  normal  coordinate  (vertical 
coordinate).  The  axes  are  scaled  by  wall  units. 

3.5  Theoretical  Correlations  for  Wall-Bounded  Turbulence 

The  modeling  approaches  described  in  section  3.4  are  clearly  not  applicable  to  wall-bounded  turbu¬ 
lence.  The  turbulence  cannot  be  considered  isotropic  or  homogeneous,  and  the  assumption  of  filter 
widths  in  the  inertial  range  are  clearly  not  justified.  However,  we  do  have  a  number  of  tools  avail¬ 
able  for  modeling  the  correlations  in  wall  bounded  flows.  In  particular,  there  is  the  representation 
theory  developed  by  Procaccia’s  group  [4,  3,  1 1]  to  represent  anisotropy  and  the  similarity  theory 
of  Oberlack  for  the  log-layer.  Further  the  quasi-normal  approximation  may  still  be  valid,  and  the 
Navier-Stokes  equations  constrain  the  possible  properties  of  the  two-point  third-order  correlation. 
Finally,  it  is  still  the  case  the  I2, 1 3  and  75  are  correlations  between  the  LES  state  variables,  so  they 
are  candidates  to  be  determined  dynamically. 

To  date,  the  application  of  the  theory  of  Procaccia  and  of  Oberlack  as  well  as  the  quasi-normal  ap¬ 
proximation  have  been  explored  by  appeal  to  DNS  channel  data.  The  results  of  these  investigations 
are  given  in  the  following  subsections. 


3.5.1  Applicability  of  the  Quasi-Normal  Approximation 


To  determine  whether  the  quasi-normal  approximation  is  applicable  in  the  near-wall  of  a  wall- 
bounded  flow,  it  was  directly  tested  in  a  turbulent  channel  flow  using  DNS  data  at  ReT  =  590. 
Various  restrictions  of  the  full  four-point  fourth  ranked  tensor  were  computed  directly  from  the 
data  and  via  the  quasi-normal  approximation  and  the  second  order  correlations  computed  from 
the  DNS.  The  results  suggested  that  the  approximation  is  remarkably  accurate  except  for  a  region 
very  near  the  wall  ( y+  <  50).  An  example  result  is  shown  in  figure  7.  At  high  Re,  the  region 
with  y+  <  50  will  need  to  be  treated  as  part  of  the  boundary  treatment  discussed  in  section  3.3, 
so  the  quasi-normal  approximation  appears  to  be  valid  where  the  interior  models  will  need  to  be 
formulated.  Details  of  the  test  of  the  quasi-normal  approximation  are  provided  in  [14],  a  preprint 
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of  which  appears  as  Appendix  B. 


3.5.2  Procaccia  and  Oberlack  Theory 

The  main  objective  of  the  study  is  to  get  a  representation  for  the  difference  between  the  two 
point  velocity  fluctuation  correlation  and  the  Reynolds  stress  ARap  (x,  r)  =  (ua(x)up(-x.  +  r))  — 
(ua(x)up(x.))  in  the  different  regions  of  the  turbulent  channel  flow.  The  particular  representation 
to  be  looked  at  is  given  by  [4],  and  has  the  following  form 

OO  l  Qmax{l) 

ARap  (x,  r )  =  ]£  E  S  c</«™  (x> r)  Btp  (?)  (23) 

l—o  m—-l  q—l 

Here,  the  basis  tensors  Bqlm  form  an  irreducible  representation  of  the  rotation  group,  implying 
that  applying  a  rotational  transformation  to  a  basis  tensor  gives  back  a  linear  combination  of  basis 
tensors  with  the  same  /  indices.  In  the  special  case  of  homogeneous  turbulence,  it  has  been  shown 
by  [4]  that  this  quality  of  the  basis  tensors  along  with  the  isotropy  of  the  governing  equations  of 
the  two  point  correlations  implies  that  it  is  possible  for  the  coefficients  cqim  (x,  r)  having  the  same 
l  to  be  proportional  to  Based  on  experimental  and  theoretical  evidence  (e.g.  [4],  [11] )  that 
£(l)  increases  with  l,  it  can  be  justified  that  Bqlm  with  higher  l  should  have  lesser  importance  as 
the  separation  r  goes  down.  This  property  is  obviously  useful,  because  it  gives  a  clear  hierarchy  of 
basis  tensors  for  the  representation  of  any  two-point  correlation  tensor. 

In  the  log-layer  of  the  channel,  [13]  has  shown  that  Rap  (x,  r)  depends  on  the  similarity  variable 
r/y,  where  y  is  the  distance  from  the  channel  wall.  Thus,  the  two-point  correlation  at  only  one  y 
in  the  log-layer  is  required  to  get  the  representation  in  the  whole  log-layer. 

Some  of  the  objectives  of  this  study  are  to  calculate  the  components  cgim  (x,  r)  of  the  basis  tensors 
and 

•  Verify  whether  it  decreases  in  importance  with  l 

•  Find  out  the  number  of  I’s  up  to  which  it  needs  to  be  retained 

•  Check  if  in  the  log-layer  they  depend  only  on  r /y  and  if  this  self-similarity  can  be  used  to  sim¬ 
plify  the  representation 


Properties  of  the  full  representation:  The  representation  of  the  two  point  correlation  containing 
terms  only  up  to  l  =  Imox  is  given  by 

Imax  l  Qmax{l) 

ARlpax)  (x,  r)  =  jr  X)  E  (x,  r)  Bq$  (r )  (24) 

1=0  7n=—/  q=  1 

The  fractional  error  between  AR)*ma:c)  and  the  exact  AR  is  given  by 


fymax  (AR^moil  —  AR|AR(,mai)  —  AR)  r2dr 
JomaT  (AR|AR)  r7dr 


(25) 


Variation  of  anisotropy  with  distance  from  wall 


Figure  8:  Relative  error  in  the  representation  of  A R  (left)  as  a  function  of  truncation  index  l  for 
various  y+,  and  (right)  when  truncated  at  l  =  0  (isotropic  term)  as  a  function  of  y+.  The  error  is 
measured  in  DNS  of  turbulent  channel  flow  at  Rer  =  940. 

With  Tfnax/h  varying  from  0.05  near  the  wall  to  0.2  at  the  center.  The  inner  product  between  any 
two  tensors  G  and  H  is  given  by 


2tt  7 r 


(G  (r)  |  H  (r)>  =  J  J  G\p  (r,  6 , <f>)  Haf3  (r,  0,  <f>)  sin  dd9d<t> 


o  o 


(26) 


The  two-point  correlations  have  been  obtained  for  channel  flow  of  ReT  =  934.  Figure  8(a)  com¬ 
pares  the  values  of  Elmax  for  various  y+  ranging  from  near  wall  to  the  center  of  the  channel.  Clearly, 
for  all  points,  a  very  good  representation  can  be  obtained  for  values  of  lmax  =  4  while  for  lmax  =  2 
we  get  a  maximum  error  of  close  to  20%.  The  odd  l  values  do  not  contribute  as  much  as  the  even  l 
values.  Therefore  if  we  include  only  the  l  =  0  and  l  =  2  components,  then  for  homogeneous  turbu¬ 
lence  we  get  10  independent  non-zero  cqim  components,  while  for  inhomogeneous  turbulence  we 
get  19  independent  non-zero  cqim  components.  We  can  also  see  that  the  curves  for  90  <  y+  <  180 
collapse,  due  to  the  fact  that  the  correlation  is  self-similar  in  the  log-layer.  Figure  8(b)  shows  the 
value  of  E°  for  different  y+  values,  i.e.  the  percentage  error  in  the  correlation  if  only  the  isotropic 
part  of  the  correlation  were  used  for  it’s  representation.  Clearly,  even  for  y+  values  close  to  the 
wall  the  isotropic  part  represents  the  majority  of  the  correlation.  Again,  evidence  of  self-similarity 
can  be  seen  by  flattening  of  the  curve  between  90  <  y+  <  180. 

Self-similarity  in  the  log-layer: 

Now  we  try  to  verify  the  self-similarity  hypothesis  in  the  log-layer.  Clearly,  if  Rap{y,  r)  depends 
only  on  r/y  then  cqim(y,  r)  should  depend  on  r/y.  Figure  9(a)  shows  the  curve  c42o(y,  r)  plotted 
against  r/y  for  different  y  values  in  the  log-layer.  We  can  clearly  see  very  good  self-similarity  in 
this  case.  However,  this  is  not  seen  in  all  the  components,  and  a  significant  proportion  of  them 
do  not  collapse  on  a  single  curve.  An  example  of  this  is  seen  in  figure  9(b)  for  cno(y,  r).  This 
implies  that  the  correlation  can  be  separated  into  self-similar  and  non  self-similar  parts.  It  is  not 
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1  =  2m  =  Oq  =  4 


1=  2m-  0  q  =  1 


Figure  9:  Angular  expansion  coefficients  cqim  for  the  indicated  values  of  the  indices.  Coefficients 
are  plotted  as  functions  of  r/y,  consistent  with  similarity  scaling  of  Oberlack  [13]. 


clear  whether  there  is  a  unique  decomposition  which  can  actually  separate  these  parts,  since  it  is 
not  possible  to  define  a  unique  non  self-similar  component.  Instead,  we  look  at  some  way  to  reduce 
the  size  of  the  representation,  such  that  for  each  ( l ,  m)  subspace  we  choose  an  optimal  tensor  and  a 
self  similar  scalar  component  which  minimizes  the  global  error  between  the  reduced  representation 
and  the  two-point  correlation  across  the  log-layer.  Thus,  we  look  for  a  representation  of  the  form 

AiC'  (n)  =  £  a,”‘(v)s§(n)  (27) 

lm 

where  r?  =  r/y  and 

s'S(n)  =  Y,}q,mBi?(n)  m 

9 

Here  fqlm  is  a  set  of  scalar  constants  and  is  one  of  the  variables  which  need  to  be  optimized,  the 
other  being  alm(rj).  To  formulate  the  optimization,  we  first  get  the  orthonormal  basis  of  Bqlm 
within  the  (/,  m)  subspace  using  the  gram-schmidt  orthonormalization  procedure.  Therefore,  we 
get 

B,qlm(fi)  =  (29) 

9' 

so  that  (B,qlm  |  B,q'lm^  =  Sqq> .  The  exact  representation  is  now  given  by 


ARa0(x,r)  =  £c'’,m(x,r 

qlm 

(30) 

We  can  redefine  the  optimal  tensors  S(m  as 

9 

(31) 

with  the  constraint 

Y^9glmgqlm  =  1 

(32) 

9 


15 


in  order  to  make  sure  that  only  alm{rj)  contains  information  about  the  magnitude  of  the  component. 
We  need  to  minimize  the  error  between  AR(x,r)  and  AR^^)  over  the  range  of  the  log  layer 
yo  <  y  <  Vi-  This  error  is  given  by 


1  Vl  t 

F  =  ——  1 1  (AR (y,  yrj)  -  AR°pt{rj)  |  AR (y,  yrj)  -  AR^fa))  rj2drjdy  (33) 

Vl  V°yo  0 


Using  the  expansions  given  by  equations  (23)  and  (27)  the  inner  product  inside  the  integral  is 
given  by  a  known  functional  of  the  coefficients  {alm(j])},  {gqlm}  and  {cqlm(y,  yrj)}.  The  problem 
is  solved  using  the  method  of  Lagrange  multipliers,  by  minimizing 

Im  q 


This  reduces  to  solving  for 


d 

dgqlm 


lm  q 


=  0 


(34) 


d 

dalm 


F  +  E'WEs,'m«‘"’n-i) 

lm  q 


=  0 


The  sets  of  {gqlm}  which  satisfy  the  above  equations  are  given  by  the  eigenvalue  problem 

EW‘*  =  Wm 

qf 

Here 

i 

Crf,m  =  / 

0 


(35) 


(36) 


(37) 


and 

i  yi 

Iqlm(y)  = - /  c'qlm(y,  yrj)dy  (38) 

yi-yoJo 

Where  {c,qlm(y,  r)}  can  be  obtained  from  {cqlm(y,r)}  using  (23),  (29)  and  (30).  Thus  {gqlm}  is 
given  by  the  eigenvectors  of  (36),  with  the  most  optimal  set  given  by  the  eigenvector  corresponding 
to  the  maximum  A*m.  The  most  optimal  set  of  {<rm(7/)}  is  then  given  by 

alTn(y)  =  '£iqlm(v)9glm  (39) 


Table  (1)  shows  the  degree  to  which  the  self-similar  representation  matches  the  two-point  correla¬ 
tion  in  the  log-layer.  Fi  shows  the  fraction  of  the  correlation  represented  by  the  individual  mode 
and  F2  is  the  fraction  of  correlation  represented  by  the  sum  of  the  first  i  modes.  Clearly,  the  self¬ 
similar  representation  is  a  very  good  approximation  of  the  correlation  across  the  log-layer.  The 
isotropic  tensor  itself  forms  the  largest  part  of  the  correlation,  while  the  l  =  2  components  come  in 
next.  Therefore,  a  self-similar  representation  using  just  the  l  =  0  and  1  =  2  components  should  be 
good  enough,  even  though  the  region  is  highly  inhomogeneous.  Figure  10  shows  the  comparison 
of  the  fully  reconstructed  self-similar  representation  with  correlations  in  a  higher  Reynolds  number 
flow  taken  from  experiments  [6].  The  self-similar  representation  gives  a  reasonable  estimate  of  the 
streamwise  correlation  across  the  log-layer. 
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0.2 


Figure  10:  Comparison  between  the  self-similar  part  of  Ru{y, r)  calculated  from  DNS  with  actual 
correlations  in  channel  flow  with  ReT  =  2433  (Experimental  data),  (a)  shows  the  self-similar 
part  of  the  correlation  calculated  from  DNS,  (b)  &  (c)  show  correlation  taken  from  experiment  for 
y/h  =  0.15  and  y/h  =  0.25,  both  in  the  log-layer.  The  correlations  are  in  the  x-y  plane,  and  the 
separations  in  the  x  and  y  directions  have  been  normalized  with  respect  to  distance  from  the  wall 


Table  1:  The  modes  in  the  self-similar  representation  sorted  according  to  decreasing  order  of  im- 

II  £  o,<i')m<<'>(?})Si(<,>n,(<'> (t})|| 

portance.  F,  =  and  IAH<Mi, - ,  where  ||M(y<r?)|P  - 


/  /  (M|M)  rfir/dy 

yo  0 


l|AR(y,yr;)|| 


1  l(i)  m(i) 

1  0  0 

2  2  1 

3  2  0 

4  2  2 

5  3  1 

6  4  0 

7  1  1 

8  1  0 

9  4  3 

10  3  2 


Fx 

f2 

0.815 

0.903 

0.066 

0.939 

0.045 

0.963 

0.024 

0.975 

0.014 

0.982 

0.012 

0.988 

0.007 

0.992 

0.003 

0.993 

0.002 

0.994 

0.001 

0.995 
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Abstract 

The  feasibility  of  an  optimal  finite-volume  large-eddy  simulation  (LES)  model  for  isotropic  tur¬ 
bulence  is  evaluated.  This  modeling  approach  is  based  on  the  approximation  of  the  ideal  LES 
by  a  stochastic  estimate  of  the  fluxes  in  a  finite-volume  representation  of  the  Navier-Stokes  equa¬ 
tion.  Stochastic  estimation  of  the  fluxes  allows  for  the  simultaneous  treatment  of  Navier-Stokes, 
discretization  and  subgrid  effects,  yielding  a  compact,  yet  accurate  scheme  for  the  large  eddy  sim¬ 
ulation  of  isotropic  turbulence.  Both  global  and  local  models  based  on  optimal  finite-volume  LES 
are  developed  and  used  in  a  priori  tests  guiding  the  choice  of  stencil  geometry  and  model  inputs. 
The  most  promising  models  in  the  a  priori  exams  are  tested  in  actual  simulations  (i.e.  a  posteri¬ 
ori)  and  the  results  compared  with  those  for  filtered  DNS  and  the  dynamic  Smagor insky  model. 
The  a  posteriori  performance  of  the  optimal  finite-volume  LES  models,  evaluated  by  the  energy 
spectrum  and  third-order  structure  function,  is  superior  to  that  of  the  dynamic  Smagorinsky  model 
on  a  coarse  grid.  While  applicability  to  other  cases  is  currently  limited  by  the  dependence  of  the 
present  approach  on  direct  numerical  simulation  (DNS)  statistical  data,  research  is  underway  to 
remove  this  requirement. 
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I.  INTRODUCTION 


In  large-eddy  simulation  (LES),  one  simulates  turbulence  by  resolving  only  the  large  scales 
while  modeling  the  small  ones.  This  approach  is  justified  by  the  observation  that  the  large 
scales  usually  dominate  the  dynamics  of  the  flow,  whereas  the  small  scales  are  important 
only  to  the  extent  to  which  they  affect  the  large  scales.  The  small  scales  of  turbulence  are 
believed  to  be  more  universal  and,  therefore,  more  amenable  to  modeling,  which  should 
allow  for  a  prediction  method  that  can  be  used,  with  little  change,  in  a  wide  range  of  flows. 
For  reviews  of  LES  see  Rogallo  &  Moin  [1] ,  Lesieur  &  Metais  [2]  and  Meneveau  &  Katz  [3] . 
Despite  its  great  promise  as  a  turbulence  prediction  technique,  current  LES  practice  suffers 
from  many  shortcomings  (e.g.  near-wall  modeling,  and  effect  of  numerical  discretization 
errors)  which  hamper  the  use  of  LES  as  a  reliable  engineering  prediction  tool.  Several  active 
research  efforts  [4-9]  are  underway  to  improve  LES  modeling,  and,  in  particular,  optimal 
LES  models  [10,  11]  as  those  reported  here  are  being  pursued  to  address  these  issues. 

Development  of  optimal  LES  has  benefited  from  a  careful  reexamination  of  the  formalism 
underlying  LES.  Because  of  the  continuous  range  of  scales  present  in  a  high  Reynolds  number 
turbulent  flow,  it  is  necessary  to  define  with  some  precision  the  large  scales  to  be  simulated 
and  the  small  scales  to  be  discarded  in  an  LES.  This  is  accomplished  through  spatial  filtering. 
Filtered  quantities,  denoted  by  (.),  represent  the  large  (or  resolved)  scales,  and  can  be  defined 
through  the  operation  of  a  filter  kernel  on  the  original  variable.  As  an  example,  the  filtered 
velocity  field  Ui  is  defined  through  the  filter  kernel  g  (in  general  g  could  be  a  distribution) 
operating  on  the  velocity  up. 

Ui(  x)  =  J  g(x',x)ui(x')dx'  (1) 


The  small  scales  in  the  LES,  e.g.  u[  =  Ui  -  uiy  are  commonly  called  subgrid  scales. 

The  standard  LES  practice  is  to  apply  the  filter  to  the  Navier-Stokes  and  continuity 
equations,  yielding  the  governing  equations  for  the  resolved  scales: 


dui 

~dt 


dui 

dxi 


+  C  =  0 


dujUj  dp  |  1  dV  M 

dxj  dxi  R edxjdxj  1 


(2) 

(3) 


2 


where 


(4) 

(5) 


Mi  =  -pi-  +  Ci 

OXj 

Tij  =  U{Uj  -  U{Uj 


and  C  and  Ci  are  commutation  terms  that  arise  when  the  filter  does  not  commute  with 
spatial  differentiation.  The  term  is  the  subgrid  stress  term,  which  represents  the  effect 
of  the  subgrid  scales  on  the  resolved  scales. 

Most  current  LES  formulations  use  discretized  versions  of  (2)  and  (3)  as  evolution  equa¬ 
tions,  with  the  modeling  effort  focused  on  determining  the  correct  representation  of  the 
subgrid  stress  term  r^.  These  subgrid  stress  models,  such  as  the  Smagorinsky  [12],  dynamic 
[13],  scale-similarity  [6,  14],  stretched  vortex  [5]  and  defiltering  models  [8,  15-17]  have  been 
used  to  simulate  a  variety  of  turbulent  flows,  but  are  derived  assuming  that  errors  due  to 
numerical  discretization  of  the  remaining  terms  in  the  evolution  equation  are  negligible. 

The  introduction  of  numerical  discretization  in  LES  raises  both  a  fundamental  and  practi¬ 
cal  issue.  The  fundamental  issue  arises  because  continuous  filters  used  in  LES  are  commonly 
invertible,  or  nearly  so.  For  example,  the  Gaussian  filter  is  formally  invertible,  and  a  top-hat 
filter  can  be  inverted  with  the  addition  of  only  boundary  information.  With  invertible  filters 
and  in  the  absence  of  numerical  discretization,  the  filtered  velocity  can  be  evolved  exactly  by 
inverse  filtering,  computing  the  terms  in  the  Navier-Stokes  equations,  then  filtering  the  re¬ 
sult.  However,  numerical  discretization  actually  eliminates  critical  small-scale  information, 
making  this  impossible,  and  thus  it  is  the  presence  of  numerical  discretization  that  makes 
LES  modeling  necessary  in  these  cases.  Because  of  this,  it  seems  desirable  that  LES  models 
be  formulated  in  the  context  of  the  numerical  representation  with  which  they  will  be  used. 

As  a  practical  matter,  the  numerical  discretization  errors  introduced  by  standard  low- 
order  numerical  schemes  are  of  the  same  order  as  the  LES  model  terms  [18-20],  leading  to 
solutions  that  are  contaminated  by  numerical  errors.  Two  different  approaches  are  currently 
used  to  avoid  such  large  numerical  errors:  high-order  discretization  schemes  (or  spectral 
methods)  and/or  grid  refinement  while  keeping  the  filter  width  A  constant.  Both  of  these 
approaches  have  been  successfully  applied  to  yield  accurate  LES  solutions  [21,  22],  but 
their  use  is  limited  by  practical  considerations.  The  application  of  high-order  methods  is 
usually  difficult  in  complex  geometries,  while  the  cost  of  adequate  grid  refinement  is  generally 
prohibitive.  Since  neither  of  these  options  is  acceptable  in  complex  turbulent  flows,  another 
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approach  seems  necessary.  In  particular,  if  subgrid  models  could  be  formulated  to  account 
for  the  numerical  discretization,  accurate  low-cost  LES  would  be  possible. 

In  the  development  of  optimal  LES,  the  above  considerations  have  led  to  an  approach 
in  which  the  numerical  discretization  is  considered  to  be  part  of  a  generalized  “filter.”  In 
this  way,  the  filter  is  viewed  as  a  mapping  from  the  infinite  dimensional  space  of  Navier- 
Stokes  solutions  to  a  finite  dimensional  LES  space,  which  can  be  represented  computationally 
without  further  truncation.  The  many-to-one  nature  of  the  mapping  leads  to  an  analysis 
of  the  subgrid  turbulence  and  its  effects  as  stochastic,  and  the  use  of  stochastic  modeling 
tools  in  the  development  of  optimal  LES  models,  which  are  deterministic  [10].  This  filtering 
approach  also  leads  to  other  differences  in  the  current  development  relative  to  other  LES 
formulations.  For  example,  in  our  formulation,  there  is  no  distinction  between  the  grid 
scale  and  the  filter  scale,  or  between  numerical  errors  and  modeling  errors.  It  would  also 
be  possible  to  formulate  optimal  LES  based  on  appropriate  non-invertible  filters  which  are 
distinct  from  the  numerical  discretization,  but  that  is  not  the  approach  taken  here. 

In  this  paper,  we  examine  the  feasibility  of  using  optimal  LES  techniques  to  develop 
models  for  use  with  a  coarse  finite- volume  representation  of  turbulence.  Finite  volume  rep¬ 
resentations  are  used  because  of  their  general  applicability  to  flows  with  complex  geometries, 
in  contrast  to  the  spectral  representations  used  in  previous  optimal  LES  studies  [10,  11].  To 
this  end,  an  optimal  finite-volume  LES  model  that  simultaneously  represents  Navier-Stokes 
dynamics,  numerical  and  subgrid  effects  is  constructed  and  evaluated.  It  is  important  to 
point  out  that,  in  the  current  study,  the  correlations  needed  for  the  optimal  modeling  ap¬ 
proach  are  supplied  from  DNS  statistical  data,  rendering  the  model  specific  to  the  isotropic 
flow  for  which  it  is  constructed.  This  is  not  an  inherent  shortcoming  of  the  optimal  LES 
approach,  but  rather  a  choice  which  was  made  to  evaluate  the  procedure  in  isolation  of 
further  modeling  which  is  necessary  for  a  complete  and  general  LES  model.  While  research 
is  underway  to  generalize  the  optimal  finite-volume  LES  model,  specifically  the  modeling  of 
the  correlations  which  are  currently  computed  from  DNS  data,  it  is  important  to  verify  if 
the  optimal  procedure,  by  itself,  is  adequate  for  constructing  the  compact,  accurate  models 
proposed  above.  A  promising  approach  to  modeling  the  required  correlations  is  discussed 
briefly  in  section  IV. 
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II.  OPTIMAL  FINITE- VOLUME  LES 


The  existence  of  a  model  defining  the  limiting  accuracy  of  an  LES  simulation  has  been 
noted  previously  [10,  23].  This  ideal  LES  model  reproduces  all  single-time,  multi-point 
statistics  exactly  and  minimizes  the  error  of  the  large-scale  dynamics.  The  evolution  of  this 
ideal  LES  is  defined  by  the  following  conditional  average: 


dio 
d  t 


(6) 


where  w  represents  an  LES  field  and  u  represents  a  real  turbulent  field  (notice  that  these 
are  entire  fields).  This  expression  for  the  ideal  LES  is  perfectly  general.  In  the  case  when 
the  filter  is  invertible,  ^  is  not  stochastic,  and  the  ideal  LES  evolution  reduces  to  ^ 
indicating  that  the  filtered  field  evolution  can  be  determined  exactly  in  terms  of  just  the 
filtered  field.  Such  a  simulation  is  equivalent  to  a  DNS,  with  the  same  resolution  require¬ 
ments.  Only  when  the  filter  is  a  many-to-one  mapping  as  described  in  section  I  does  this 
formalism  correspond  to  an  LES,  and  this  is  the  case  of  interest  here. 

Unfortunately,  direct  computation  of  the  ideal  model  is  not  practical,  because  of  the 
immense  amount  of  statistical  information  contained  in  (6).  Optimal  LES  is  a  formal  ap¬ 
proximation  of  the  ideal  LES  using  stochastic  estimation,  a  well-known  technique  for  ap¬ 
proximating  the  conditional  average  with  a  manageable  amount  of  data  [24-26].  Optimal 
LES  has  been  previously  applied  to  model  the  subgrid  term  with  spectral  numerical  repre¬ 
sentations  of  the  Navier-Stokes-like  terms  in  the  filtered  equation[10,  11].  In  these  previous 
studies,  the  focus  was  the  analysis  of  the  resulting  optimal  subgrid  model  in  order  to  identify 
the  fundamental  properties  a  subgrid  model  should  have  in  isotropic  [10]  and  wall-bounded 
[11]  turbulence.  The  authors  made  no  attempt  to  construct  a  general  optimal  LES  model 
for  use  in  a  range  of  complex  turbulent  flows. 

One  of  the  primary  impediments  to  the  broad  application  of  the  optimal  LES  approach 
developed  in  these  previous  studies  is  their  reliance  on  spectral  numerical  representations, 
which  are  not  easily  applied  to  flows  in  complex  geometries.  In  the  present  study,  the 
applicability  of  optimal  LES  to  finite  volume  representations  is  explored.  This  is  more 
challenging  than  spectral  methods,  since  exact  derivative  operators  cannot  be  defined.  To 
avoid  the  numerical  errors  inherent  in  a  finite-volume  representation  of  the  Navier-Stokes 
terms,  the  optimal  modeling  procedure  is  applied  to  the  filtered  Navier-Stokes  equation  as 


5 


a  whole,  rather  than  just  the  subgrid  stress  as  was  done  in  [10,  11].  A  generally  applicable 
optimal  finite-volume  model  requires  modeling  of  the  small-separation  second-,  third-  and 
fourth-order  velocity  correlations  appearing  in  section  II B  3.  For  the  current  study,  this 
requirement  is  being  bypassed  through  the  use  of  DNS  statistical  data,  allowing  us  to  focus 
on  the  viability  of  the  finite- volume  optimal  LES  modeling  procedure. 


A.  Finite- volume  formulation  of  LES 


To  evaluate  the  feasibility  of  applying  the  optimal  LES  procedure  directly  to  a  finite 
volume  representation,  a  coarse  finite- volume  representation  of  homogeneous  isotropic  tur¬ 
bulence  is  explored.  In  this  case,  the  domain  is  a  periodic  cube  of  dimension  2tt,  with  the 
numerical  grid  for  each  velocity  component  composed  of  N 3  cubes  of  dimension  A  =  2-n/N. 
To  allow  for  staggered  as  well  as  collocated  grids,  the  cubes  are  labeled  as  Ctji,  each  with 
center  located  at  xji  =  xj0  +  iA,  where  xj0  is  the  origin  of  each  grid.  For  collocated  grids 
xj0  =  0,  and  for  staggered  grids  xj0  =  -f  e5',  where  ej  is  the  unit  vector  in  the  j  direction. 
The  j  indices  refer  to  the  grid  associated  with  the  jth  component  of  velocity,  while  i  is  a 
vector  of  integers  which  establishes  the  position  of  the  cube. 

The  finite- volume  filter  kernel,  g  (x  -  £)  as  defined  in  (1)  is  a  top-hat  filter  sampled  on  a 
grid  and  is  given  by: 

gji  (x  -  0  =  S(xji  -  x)G  (*!  -  &)  G  (x2  -  6)  G  (a*  -  6)  (7) 


where  8  is  the  Dirac  delta  function  and  G  (xi  -  &)  is  the  one-dimensional  box  or  top-hat 
filter: 


<?(*<-&)  = 


i  if  |*<-{<|  <4, 


(8) 


y  0  otherwise. 

Note  that,  to  support  staggered  grids,  the  filter  kernel  depends  on  which  velocity  component 
is  being  filtered. 

The  finite-volume  filtered  velocity  is,  therefore,  defined  as: 


wj  =  ^Jn  ui(x)dx •  (9) 

The  evolution  equations  for  the  filtered  velocity  are  found  by  filtering  the  Navier-Stokes 
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equations  with  g,  yielding  the  control-volume  form  of  these  equations: 


The  divergence  theorem  is  used  to  convert  the  volume  integrals  to  surface  integrals: 


(10) 

(11) 


A3ir + 1,  u^dx’ = -  i,  + h  L„  Untdx'  (12) 

/  UjTijdx'  =  0  (13) 

J  d£lji 

where  rij  is  the  outward  pointing  normal  of  the  volume  boundary  dttji. 

The  evolution  of  the  LES  field  now  depends  only  on  the  surface  fluxes  in  each  volume. 
In  a  standard  finite- volume  scheme,  the  fluxes  would  be  approximated  numerically  through 
reconstruction  techniques,  yielding  approximations  which  converge  as  the  grid  size  tends  to 
zero,  A  — *  0.  However,  such  convergence  considerations  are  not  applicable,  for  the  LES 
considered  here,  since  A  >  ?/  (r|  is  the  Kolmogorov  scale).  Instead,  the  convective  and 
diffusive  fluxes  in  equation  (12)  will  be  modeled  using  stochastic  estimation.  One  could 
also  subject  the  pressure  flux  term  to  the  same  modeling  technique,  but,  as  pointed  out 
by  Langford  &  Moser  [27],  such  a  representation  would  be  inherently  non-local  and  thus 
impractical  in  complex  applications.  An  alternative  is  to  use  an  approximate  continuity 
constraint  to  determine  the  pressure  fluctuation,  as  in  traditional  finite-volume  methods. 
Langford  &  Moser  [27]  have  shown  that  an  optimal  scheme  for  pressure  is  only  slightly  more 
accurate  than  the  traditional  schemes,  especially  in  the  staggered  grid  case. 

The  task  of  modeling  the  surface  fluxes  appearing  in  the  exact  LES  equations  (12)  can 
be  simplified  considerably.  The  flux  across  a  particular  face  appears  in  the  momentum 
equation  of  two  neighboring  cells;  to  ensure  conservation,  the  two  cells  must  have  identical 
representations  of  the  flux.  In  isotropic  turbulence,  homogeneity  dictates  that  a  model  for 
the  flux  through  a  given  face  must  not  depend  on  the  location  of  the  face.  These  two 
requirements  imply  that  a  model  need  only  be  provided  for  three  of  the  six  faces  of  the 
cube. 

Further  simplifications  are  introduced  by  isotropy.  The  symmetry  of  the  grid  requires 
that  the  equations  for  the  LES  evolution  be  invariant  to  reflections  and  rotations  in  the 
coordinate  system.  This  basic  requirement  imposes  a  symmetry  on  the  flux  models,  and 
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implies  that  each  flux  term  can  be  constructed  from  rotations  and  reflections  of  four  (4) 
basic  fluxes.  The  terms  M^w,  M*w,  Mfjw/dz  and  M£u/dz  are  used  to  denote  these  basic 
fluxes,  which  are  averages  of  ww,  uw,  dw/dz,  and  du/dz  (for  the  remainder  of  this  paper, 
u,  v,  and  w  are  used  to  denote  components  of  w),  respectively,  over  a  face  of  dimension 
A  x  A  oriented  with  normal  in  the  ^-direction  and  located  at  the  top  (positive  z  side)  of 
the  ith  volume  element.  For  Mww  and  Mdw/dz,  the  relevant  volume  element  is  for  the  w- 
component  of  velocity,  and  for  the  other  fluxes,  the  relevant  volume  is  for  the  u-component 
of  velocity.  The  four  fluxes  to  be  modeled  are 


1 

r4l+ a/2 

/■xf+A/2 

(14) 

AT  = 

WW 

X 

A  o 

/ 

f 

uiu-z&x&y 

lxf-A/2  J 

xf-A/2 

1 

/•x“+A/2 

/•x«+A/2 

(15) 

Ml  = 

uw 

X 

AO 

/ 

f 

u3u3dxdy 

A2  v 

lxtf-A/2  J 

Xj*— A/2 

1 

/»X2*+A/2 

/■xf+A/2 

(16) 

Mdw/dz  z 

X 

=  A*, 

Jxf-A/2  v 

/xf-A/2 

^dxdy 

oz 

Mdu/dz  - 

1 

"A*. 

/.x»+A/2 

Jig-  A/2  . 

/■xf+A/2 

/xf-A/2 

dii3A  , 
dxdy 
oz 

(17) 

To  actually  perform  a  simulation,  the  models  for  these  fluxes  must  be  rotated  and  reflected 
appropriately  to  generate  all  the  terms  in  (12). 

B.  Stochastic  estimation  technique 

In  stochastic  estimation,  which  is  the  basis  for  optimal  LES  models,  one  begins  by  defining 
a  restricted  functional  form  (e.g.  linear  or  quadratic)  for  the  model.  The  stochastic  estimate 
is  then  determined  by  minimizing  the  mean-square  difference  between  the  model  and  the 
exact  quantity  being  modeled.  The  application  of  this  technique  to  the  finite  volume  fluxes 
is  discussed  below,  followed  by  a  brief  discussion  of  the  sense  in  which  such  model  can  be 
considered  to  “converge.” 
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1.  Model  definition 


The  stochastic  estimation  model  rri1  of  a  particular  flux  term  M*,  is  expressed  as  a  sum 
of  linear  and  quadratic  products  of  filtered  velocities: 

m‘  =  A +  '£  BP w‘+m  +  E  qr<m«4+”  (18) 

m  m  n 

More  elaborate  model  forms  are  possible,  but  at  least  a  quadratic  dependence  is  needed, 
since  the  convective  fluxes  are  quadratic.  Note  that  m  and  n  are  vectors  representing  the 
separation  between  the  face  being  modeled  (face  i)  and  the  velocity  term  used  in  the  sum. 
One  can  extend  the  summation  over  all  volumes,  generating  a  global  model,  or  restrict  m 
and  n  to  a  certain  stencil  of  volumes,  generating  a  local  model. 

The  geometry  of  the  stencils  for  local  models  is  described  through  the  following  notation: 

x  x  Nz  -Nvxx  Nvy  x  Nvz  -  Nfi  x  N™  x  Nfi 

where  Nx,  Ny  and  Nz  define  the  number  of  volumes  in  the  stencil  in  each  direction.  In  our 
convention,  the  fluxes  are  estimated  on  a  face  with  surface  normal  pointing  in  the  positive 
^-direction,  with  the  stencil  centered  about  this  face.  The  u,  v  and  w  superscripts  specify 
the  component  of  velocity  for  which  the  stencil  is  used.  This  distinction  is  necessary  since, 
in  general,  different  velocity  components  will  have  different  stencils  for  a  given  flux.  When 
referring  to  a  stencil  for  a  single  velocity  component,  the  simpler  notation 

Nxx  Nyx  Nz 

can  be  used.  In  Fig.  1,  examples  of  single  velocity  stencils  are  shown. 

2.  Defining  the  minimization 

In  the  stochastic  estimation  technique,  the  model  coefficients  A ,  B  and  C  are  determined 
by  minimizing  the  mean-square  error  with  respect  to  the  model’s  inputs  or  events  (u>j+Tn 
and  wj+mit4+n  in  the  model  defined  by  equation  (18)).  The  stochastic  estimate  is  intended 
to  approximate  the  conditional  average  constituting  the  ideal  LES  (in  this  case  m*deal  = 
( Mi\u  =  w)),  so  it  is  the  mean-square  error  in  reproducing  m*deal  that  is  to  be  minimized, 
even  though  mfdeal  is  not  known.  However,  it  can  be  shown  [11]  that  this  error  is  also 
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minimized  when  the  mean-square  difference  between  M*  and  mt  is  minimized.  Thus,  it  is 
this  mean-square  difference 

d1  =  M*  -  m*  (19) 


that  is  directly  minimized  in  developing  the  estimates.  We  refrain  from  labeling  d  as  an  error 
since  it  cannot,  even  in  principle,  be  made  to  approach  zero.  Its  mean  square  is  bounded 
from  below  by  the  variance  of  Ml  about  mfdea! .  Instead  d  will  be  called  the  “variation”  of 
the  model. 

However,  there  are  at  least  two  different  ways  to  define  the  variation  when  estimating  the 
fluxes.  The  single-flux  variation  measures  the  difference  between  the  model  and  the  flux  at 
a  single  face,  and  is  defined  as 

d\  =  Mi-  mi  (20) 


While  this  is  the  most  straightforward  definition  of  the  model’s  variation,  it  is  not  the  most 
relevant  to  the  dynamics  of  the  LES  evolution.  For  example,  if  the  flux  contributes  to  the 
momentum  equation  for  the  w  component  of  velocity,  as  for  Mww  and  Mqw/qz,  then  the 
model  appears  in  the  momentum  evolution  equation  as  follows: 


^  =  •  ■  •  +  m*  - 
df 


(21) 


The  contribution  of  the  model  to  the  variation  in  the  momentum  evolution  is  the  dual-flux 
variation  of  the  estimate,  defined  as 


dj |  =  (AT4  -  mi)  -  (Mi+(0’0>1)  -  mi+(0’0,1))  (22) 

In  general,  the  dual- flux  variation  d\\  is  different  from  the  single-flux  variation  d\,  because 
there  are  stochastic  portions  of  the  exact  flux  term  M‘  that  exactly  cancel  portions  of  the 
flux  term  on  the  opposing  face.  The  single-flux  variation  d\  includes  contribution  from  the 
cancellation,  and  may  thus  overpredict  variation  in  the  momentum  equation.  Whenever  the 
variation  is  plotted  or  tabulated,  d\\  is  used,  since  this  non-canceling  portion  of  the  variation 
is  most  relevant  to  the  dynamics  of  the  LES,  due  to  its  explicit  appearance  in  the  momentum 
equations. 

In  reporting  a  priori  variations,  it  is  clearly  d^  that  is  of  interest,  but  one  of  the  other 
variations  may  be  most  appropriate  to  minimize  to  determine  the  modeling  coefficients  A, 
B  and  C.  When  local  estimates  are  used,  somewhat  different  models  will  result  depending 
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on  which  variation  is  minimized,  and  the  differences  can  be  important.  For  example,  it 
is  shown  in  Appendix  A  that  minimizing  d\  results  in  a  model  that  conserves  momentum, 
and  correctly  predicts  dissipation  in  the  a  priori  sense.  Whereas,  minimizing  d||  results 
in  models  that  are  either  conservative  or  have  correct  a  priori  dissipation,  but  not  both. 
Local  estimates  were  computed  using  both  the  single-flux,  d\,  and  dual-flux,  d||,  variations, 
and  it  was  found  that  the  increase  in  d\\  due  to  minimizing  the  mean-square  of  d\  is  small, 
but  that  the  error  in  predicting  a  priori  dissipation  is  significant  when  the  mean-square  of 
d\\  is  minimized.  For  these  reasons,  the  local  models  developed  here  were  determined  by 
minimizing  d\. 


3.  Determining  Optimal  Estimates 


In  order  for  a  linear  estimate  to  minimize  d\,  the  following  conditions  should  be  met  [28]: 

<AP)  =  (m*)  (23) 

((M*  -  rrd)Et)  =  0  (24) 


where  E  is  the  event  vector,  consisting  of  all  events  used  in  the  model  equation  (18).  For 
an  optimal  linear  estimate  (C  =  0  in  (18))  the  event  vector  is  given  by: 


E 


=  (<ra) 


(25) 


while  the  event  vector  for  the  optimal  quadratic  model  is 

t+m 


E  = 


w 


Wj 
i+m  i+n 


(26) 


where  m  and  n  vary  over  all  separations  defined  by  the  stencil  for  a  specific  model. 

Therefore,  to  determine  the  optimal  linear  estimates  and  the  associated  variation,  the 
following  correlations  are  required: 


(MM) ,  (M<> ,  «*<) 


(27) 


The  i  superscripts  have  been  removed  since  the  correlations  are  averaged  over  all  volumes.  To 
compute  the  optimal  quadratic  estimate,  the  following  additional  correlations  are  required: 

,  «*«} ,  (28) 
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The  data  requirement  for  representing  the  entire  fourth-order  correlation  function  in  (28)  is 
prohibitive.  To  ease  the  requirement,  quadratic  products  in  the  model  (18)  are  restricted  to 
nearby  separations  (|m  —  n|  <  1). 

The  linear  optimal  estimate  is  determined  using  the  following  estimation  equations: 

(M)  =  A  (29) 

(Mw d  =  y,b r«»?>  <3°) 

m 

For  the  quadratic  optimal  estimate,  the  following  set  of  estimation  equations  must  be  solved: 

(M)  =  A  +  E  J2  °T  «*»?>  (31) 

m  n 

(Mw?)  =  Br  «<«r> + E  E  «« >  <32) 

m  m  n 

(Mwfw?)  =  A  (wfw?)  +  22Br  «<">?>  +  E  E ' CT  (33) 

m  m  n 

In  this  work,  the  correlations  are  computed  using  filtered  DNS  fields  from  a  well-resolved 
pseudo-spectral  2563  simulation  of  forced  isotropic  turbulence.  Numerical  details  of  the 
simulation  can  be  found  in  Rogallo  [29].  In  the  arbitrary  units  of  the  simulation  code,  the 
DNS  data  has  the  following  characteristics:  turbulent  kinetic  energy  q2/ 2  =  41.1,  dissipation 
e  =  62.9,  Taylor  microscale  A  =  0.209,  and  microscale  Reynolds  number  R\  =  164.  The 
three-dimensional  energy  spectrum  of  the  DNS  is  displayed  in  Fig.  2,  with  the  filter  width 
clearly  indicated.  This  is  the  same  DNS  used  by  Langford  k  Moser  [10]  for  computing  their 
optimal  subgrid  stress  model.  Note  that  due  to  the  relatively  low  Reynolds  number  of  the 
DNS,  the  filter  width  A  could  not  be  either  much  larger  or  much  smaller  in  this  flow  if  it  is 
to  remain  in  the  approximate  inertial  range. 

The  correlations  were  gathered  from  15  DNS  fields,  each  separated  by  approximately 
one-half  of  an  eddy-turnover  time.  To  increase  the  statistical  sample,  each  of  the  48  possible 
reflections  and  rotations  of  the  field  were  also  included  in  the  averaging.  Note  that  the 
mapping  from  a  single  DNS  field  to  a  filtered  field  is  not  unique,  because  the  volume  elements 
can  be  positioned  at  many  different  locations.  When  mapping  from  a  2563  DNS  field  to  a 
323  LES  field,  83  linearly  independent  choices  exist  for  the  mapping.  By  averaging  over  all 
the  possibilities,  one  ensures  the  highest  quality  correlations.  In  this  work,  a  choice  was 
made  to  consider  only  23  LES  mappings  for  each  DNS  field.  This  number  was  chosen  to 
allow  a  data  set  that  could  be  stored  and  managed  on  available  equipment. 
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4-  Convergence 


In  traditional  finite-volume  schemes  for  incompressible  flows,  the  flux  model  is  derived 
from  a  reconstruction  of  the  flux  function  from  the  cell-averaged  velocity  data.  A  reconstruc¬ 
tion  based  on  a  Taylor  series  expansion  is  usually  appropriate  because  it  can  be  assumed 
that  the  underlying  function  is  smooth  and  that  the  grid  separation,  A,  is  small  when  com¬ 
pared  to  the  smallest  lengthscales  of  the  velocity  field.  The  resulting  approximation  then 
converges  in  the  usual  sense,  in  that  the  errors  decrease  with  decreasing  A,  and  can  be  made 
arbitrarily  small,  limited  only  by  the  precision  of  the  arithmetic  used  in  the  calculations. 
For  example,  order  of  accuracy  considerations  in  traditional  finite- volume  schemes  results  in 
an  error  decaying  as  A2  (a  second  order  scheme)  for  a  1x1x2  stencil,  while  a  1x1x4  stencil 
yields  a  fourth-order  accurate  scheme. 

The  ideas  of  convergence  and  order  of  accuracy  of  the  numerical  model  must  be  reex¬ 
amined  in  the  current  optimal  finite-volume  LES  approach.  In  this  case,  A  is  by  definition 
not  small  compared  to  the  smallest  scales  of  the  velocity  field.  Furthermore,  in  the  LES 
formulation  considered  here,  the  filter  definition  and  numerical  discretization  are  identical, 
so  the  grid  size  A  cannot  be  small  compared  to  the  filtered  velocity  either.  Reconstruction 
of  the  fluxes  based  on  Taylor  series  or  other  numerical  procedure  would  clearly  be  invalid  in 
this  case,  as  are  our  usual  notions  of  convergence. 

In  optimal  finite  volume  LES,  two  senses  of  convergence  should  be  considered,  though 
neither  will  be  of  direct  utility  in  the  current  work.  First,  at  fixed  A  we  can  ask  whether 
a  sequence  of  increasingly  complex  models  converges  to  the  ideal  model.  For  the  model 
defined  in  (18),  the  complexity  can  be  increased  by  increasing  the  size  of  the  stencil.  But 
this  clearly  does  not  yield  convergence  to  the  ideal  model,  since  (18)  is  restricted  to  quadratic 
dependence  on  w.  While  this  does  not  yield  convergence  to  the  ideal,  the  way  local  models 
approach  the  corresponding  global  model  as  the  stencil  is  expanded,  is  of  interest  [30].  In 
addition,  adding  higher  order  terms,  (cubic,  fourth  order  etc.)  to  (18)  could  in  principle 
result  in  convergence  to  the  ideal,  though  this  would  be  difficult  to  prove,  and  impractical 
to  implement. 

The  second  sense  of  convergence  is  statistical.  For  a  given  non-ideal  model,  we  are 
interested  in  the  convergence  of  statistical  quantities  representative  of  the  large  scales  with 
decreasing  A.  This  is  the  sense  of  convergence  defined  by  Pope  (private  communication, 


13 


2002;  see  also  Meneveau  [31]).  In  this  case,  A  is  considered  to  always  be  much  larger  than 
the  Kolmogorov  scale  (formally  Re  — »  oo  as  A  — >  0),  so  the  usual  convergence  to  a  DNS 
does  not  apply.  To  capture  the  large  scale  statistics,  it  is  clearly  necessary  for  A  to  be 
small  compared  to  the  large  scales  of  turbulence  (which  are  of  order  q 3  / e,  where  q2  is  twice 
the  turbulent  kinetic  energy  and  e  is  the  dissipation),  thus  A e/q3  should  be  small  (it  is 
0.017  for  the  case  considered  in  section  III).  For  comparison  purposes,  the  ratio  of  A  to 
the  Kolmogorov  lengthscale,  A/77,  is  23.7,  and  the  ratio  of  A  to  the  longitudinal  integral 
lengthscale,  A /Ln,  is  0.195.  Convergence  in  this  sense  clearly  depends  on  the  statistical 
quantity  considered.  For  the  purposes  of  this  paper,  we  will  consider  the  low-wavenumber 
spectrum  and  the  large-separation  third-order  structure  function.  However,  because  of  the 
reliance  on  moderate  Reynolds  number  DNS  for  the  statistical  data  needed  in  the  optimal 
model,  it  is  not  possible  to  empirically  explore  statistical  convergence  in  the  limit  defined 
above. 

III.  RESULTS 

Diagnostics  of  LES  models  are  usually  divided  into  two  distinct  categories:  a  priori  and  a 
posteriori.  In  a  priori  tests  one  evaluates  the  model’s  performance  by  comparing  the  modeled 
terms  with  actual  terms  obtained  from  DNS  (or  experiments,  where  possible).  Notice  that 
this  kind  of  diagnostic  is  actually  an  integral  part  of  the  optimal  modeling  technique,  since 
the  optimal  model  minimizes  a  given  a  priori  variation  (in  our  case,  the  single-flux  variation 
d\  defined  by  (20)). 

A-priori  testing  of  LES  models  has  often  revealed  that  the  modeled  term  correlates  rather 
poorly  with  the  actual  term  determined  from  real  turbulence  realizations.  In  the  language 
of  the  current  study,  the  variation  associated  with  the  model  is  large.  As  was  mentioned 
in  section  II B  2,  it  is  inappropriate  to  label  this  variation  “error,”  because  much  of  this 
commonly  identified  “error”  is  actually  due  to  the  stochastic  nature  of  the  filtered  field 
evolution.  Due  to  the  many-to-one  nature  of  the  filter  mapping,  a  filtered  field  does  not 
include  enough  information  to  uniquely  determine  its  evolution,  or  for  a  finite  volume  filter, 
the  fluxes.  The  variability  in  the  possible  filtered  evolutions  is  responsible  for  the  irreducible 
variance  about  the  ideal  LES  evolution  given  by  (6). 

Ideally,  one  should  separate  the  variation  measured  a-priori  for  any  model  into  two 
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components;  the  deterministic  error  of  the  model  and  the  stochastic  variance.  The  stochastic 
variance  is  intrinsic  to  the  LES  and  is  the  variation  attained  by  the  ideal  LES  model.  This 
portion  of  the  a-priori  variation  is  only  dependent  on  the  characteristics  of  the  filter  and  the 
turbulent  flow.  The  deterministic  error  reflects  modeling  errors  and  quantifies  the  accuracy 
of  a  given  model  in  approximating  the  ideal.  Unfortunately,  these  components  of  the  a-priori 
variation  cannot  be  computed  separately.  Only  the  combined  variation  of  a  given  model  can 
be  determined,  and  it  can  give  the  false  impression  that  the  model  is  inaccurate  when  it  is 
the  stochastic  variance  that  is  large. 

In  many  models  found  in  the  literature  [32-35],  the  correlation  between  the  a  priori  and  a 
posteriori  results  has  been  poor,  since  performance  commonly  appears  to  be  unsatisfactory 
in  the  former,  but  is  acceptable  in  the  latter.  This  suggests  that  the  a-priori  results  may 
have  been  dominated  by  the  stochastic  variance  component,  which  has  no  effect  on  single¬ 
time  multi-point  statistics  [10,  23].  To  compare  LES  models  fairly  in  a-priori  evaluation, 
one  must  be  careful  to  keep  the  stochastic  portion  of  the  results  consistent.  In  the  present 
case,  the  a-priori  results  for  all  models  are  computed  with  the  same  filter  and  for  the 
same  turbulent  flow.  Therefore,  the  a  priori  variations  for  the  models  are  indicative  of 
relative  model  accuracy.  In  particular,  the  measured  difference  in  the  mean-square  variation 
associated  with  two  models  is  also  the  mean-square  difference  in  the  deterministic  error, 
which  can  be  used  as  a  guide  in  selecting  the  best  stencil  geometries  and  model  inputs. 

The  a  posteriori  tests  measure  the  performance  of  the  models  in  actual  simulations,  with 
emphasis  on  the  simulation  results,  rather  than  a  detailed  analysis  of  the  models’  accuracy. 
The  models  with  the  best  performance  in  the  a-priori  tests  will  be  evaluated  in  simulations, 
where  they  will  be  compared  with  the  dynamic  Smagorinsky  subgrid  model  and  filtered  DNS 
results. 


A.  A  Priori  Analysis  of  Global  Estimates 

Relative  root  mean-square  variation  measurements  for  the  global  estimates  are  shown  in 
Table  I.  Shown  are  the  a-priori  dual-flux  variations  dy,  which,  as  noted  before,  are  a  better 
measure  of  impact  on  the  dynamics  of  an  LES  than  the  single-flux  variations  d\.  The  RMS 
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variations  are  defined  as: 


RMS  = 


^  /(M<  -  M'+CA1))2^ 


(34) 


Note  that  the  difference  between  the  mean-square  variations  of  two  models,  8d,  is  a 
measure  of  the  difference  in  the  deterministic  modeling  error  between  these  models.  One 
can  easily  verify  the  relation  between  the  variation  and  deterministic  error  differences  by 
defining  8d  as: 


8d  =  RMS22  -  RMS2  (35) 

/(M*  -  Mi+(°’°’1'>)2^ 

where  the  1  and  2  subscripts  refer  to  two  different  models. 

Writing  the  dual-flux  variation  d\\  as  the  sum  of  the  modeling  error  e  and  stochastic 
variance  v 

d\\=e  +  v ,  (36) 


substituting  (36)  in  (35)  and  expanding  yields 

(e2)2  ~  (e2)i  +  2  (eu)2  -  2  (ev)^  +  (u2)2  -  ( v ^ 
((Mi  -  Mi+(°’0’1'))2^ 


Sd 


(37) 


Equation  (37)  can  be  simplified  by  noting  that  the  stochastic  variance  is  a  property  only  of 
the  filter  and  the  flow,  not  of  the  model,  so  that  (v2)2  =  (v2)v  Furthermore,  the  modeling 
error  and  stochastic  variance  are  orthogonal  (see  [11]);  that  is  (ev)2  =  (ev}1  =  0.  Therefore, 
8d  is  the  difference  in  mean-square  modeling  error  between  two  models 


Sd  = 


(e2)2  -  (e2), 
/(Mi  -  M^0-0-1))2^ 


(38) 


Variation  measurements  can  thus  be  used  as  a  measure  of  relative  modeling  accuracy  between 
the  optimal  models. 

These  variation  measurements  allow  one  to  make  comparisons  of  the  following:  linear  vs. 
quadratic  estimates;  staggered  vs.  collocated  grids;  local  vs.  nonlocal  quadratic  products. 
Local  quadratic  products  are  those  in  which  the  velocity  components  forming  the  quadratic 
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event  variable  are  located  in  the  same  or  overlapping  volumes.  For  nonlocal  quadratic  prod¬ 
ucts,  the  separation  between  velocity  components  may  be  greater.  One  can  also  assess  the 
importance  of  terms  which  are  not  similar  to  those  in  standard  finite- volume  schemes.  These 
terms  will  be  referred  to  as  nonstandard.  Standard  terms  are  those  that  would  arise  in  stan¬ 
dard  finite- volume  schemes:  w  velocities  for  Mqw/qz  ;  u  velocities  for  Mqu/qz  ;  w-w  velocity 
products  for  Mww\  u-w  velocity  products  for  Muw.  Estimates  having  only  standard  terms 
are  marked  with  a  f  in  Table  I.  Nonstandard  dependencies  in  the  model  might  improve  esti¬ 
mates  because  the  velocity  components  are  correlated,  so  that  velocity  components  on  which 
a  quantity  is  not  directly  dependent  (e.g.  v  in  the  the  Mww  flux)  can  provide  information 
about  that  quantity. 

The  following  conclusions  are  supported  by  the  variation  measurements: 

•  Linear  estimates  are  suitable  for  estimations  of  the  viscous  fluxes,  Mdw/dz  and  Mqu/qz  , 
but  not  for  the  convective  fluxes  (as  expected). 

•  For  the  convective  fluxes,  Mww  and  Muw  on  collocated  grids,  the  inclusion  of  nonlocal 
quadratic  products  provides  modest  (<  3%  of  the  dual  flux  term)  reduction  in  error 
when  only  standard  terms  are  used,  and  a  larger  (>  5%)  reduction  in  error  when 
nonstandard  terms  are  used. 

•  For  estimates  with  only  standard  terms,  there  is  no  clear  advantage  to  either  the 
staggered  or  collocated  grid.  When  nonstandard  terms  are  included,  the  staggered 
grid  is  better  for  the  viscous  fluxes,  and  the  collocated  grid  is  better  for  the  convective 
fluxes. 

•  For  the  viscous  fluxes,  including  nonstandard  terms  provides  a  modest  (3-4%  of  the 
dual  flux  term)  reduction  in  error.  For  the  convective  fluxes,  the  use  of  nonstandard 
terms  reduces  error  by  up  to  6%. 

For  estimates  with  only  standard  terms,  the  relative  root  mean-square  variations  are 
about  36%  for  Mdw/dz ,  51%  for  Mdu/dz,  44%  for  Mww,  and  50%  for  Muw.  These  variations 
are  large,  especially  considering  that  they  are  normalized  by  the  magnitude  of  the  total 
flux  term,  which  represents  both  the  Navier-Stokes  and  subgrid  effects.  The  large  a  priori 
variations  of  these  global  estimates  are  similar  to  those  encountered  by  Langford  &  Moser 
[10]  and  Volker,  Moser  &  Venugopal  [11]  for  the  subgrid  terms.  As  pointed  out  previously, 
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such  large  variations  are  to  be  expected  if  the  term  being  modeled  is  largely  stochastic. 
Particularly,  if  the  variance  about  the  ideal  model  dominates  the  variation,  then  a  large 
variation  will  not  impact  the  accuracy  of  an  LES. 

B.  A  Priori  Analysis  of  Local  Estimates 

While  a  global  model  will  provide  the  best  estimate  for  the  fluxes,  the  computational  cost 
associated  with  the  large  stencil  (32x32x32  volumes  for  each  velocity  component)  makes  it 
impractical  for  actual  simulations,  and  may  be  unnecessary.  Local  models  may  be  nearly  as 
accurate,  at  a  fraction  of  the  cost  of  a  global  model.  A  large  collection  of  local  models  was 
computed  for  both  staggered  and  collocated  grids.  The  complete  set  of  a  priori  results  can 
be  found  in  Langford  [30]. 

The  observations  made  in  the  case  of  global  estimates  still  hold  for  the  local  models. 
For  the  viscous  fluxes,  linear  estimates  are  sufficient.  For  convective  fluxes  the  nonlocal 
quadratic  models  result  in  the  lowest  variations.  There  is,  however,  a  dependence  on  stencil, 
flux  and  type  of  grid  on  the  effect  of  using  standard  or  nonstandard  estimates. 

The  viscous  terms  have  little  sensitivity  to  the  use  of  nonstandard  terms,  notably  in  the 
thinner  stencils  (the  width  of  the  stencil  is  related  to  the  number  of  volumes  in  the  direction 
parallel  to  the  face,  while  the  length  is  related  to  the  number  of  volumes  in  the  direction 
normal  to  the  face).  With  wider  stencils,  accuracy  improves  with  the  use  of  nonstandard 
terms,  but  the  improvement  is  small.  For  example,  the  best  linear  nonstandard  staggered 
model  for  Mdu/dz{ the  5x5x6-6x6x6-6x5x5  stencil)  shows  only  a  7%  improvement  in  the 
variation  over  the  much  thinner  linear  standard  staggered  Ixlx6-2x2x2-2xlxl  stencil. 
Note  that  the  smaller  stencil  is  composed  of  only  16  volumes,  while  the  best  model’s  stencil 
is  composed  of  516  volumes.  The  cost  of  the  model  increases  faster  than  linearly  with  the 
number  of  volumes  ( N 2,  or  Nlog(N)  if  a  Fourier  transform  is  used),  so  one  must  weigh  the 
small  accuracy  improvement  against  the  large  increase  in  cost. 

The  effect  of  nonstandard  terms  on  the  convective  fluxes  is  also  stencil  dependent.  The 
effect  on  the  thinner  stencils  is  small,  but  the  improvement  on  the  wider  stencils  is  quite 
significant,  especially  in  the  collocated  case  (up  to  19%  improvement).  The  improvement  in 
the  case  of  staggered  grids  is  not  as  important  (a  maximum  of  6%).  Once  again  one  must 
measure  the  cost  of  using  a  wider  stencil  against  the  reduction  in  error. 
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When  extending  the  stencil,  one  must  weigh  the  advantages  of  lengthening  versus  widen¬ 
ing.  It  is  clear  from  the  a  priori  results  that  lengthening  is  preferable  in  almost  every 
occasion.  In  the  viscous  flux  case,  widening  the  stencil  brought  little,  if  any,  improvement  to 
the  variation,  while  lengthening  the  stencil  from  2  to  6  volumes  decreased  the  variation  by 
as  much  as  22%.  In  the  case  of  the  convective  terms,  the  same  observation  holds,  with  the 
caveat  that  if  nonstandard  terms  are  used,  widening  of  the  stencil  does  produce  a  noticeable 
decrease  in  variation. 

The  difference  in  variation  between  the  best  local  estimates  and  the  global  estimates  is 
small.  The  variation  for  the  best  estimates  (with  stencils  which  occupy  approximately  1/25 
of  the  volume  occupied  by  the  global  stencil)  are  within  10%  of  the  global  estimates,  and 
most  of  the  improvement  in  accuracy  occurs  when  lengthening  the  stencil  to  four  (4)  volumes 
in  the  direction  normal  to  the  face.  This  indicates  that  most  of  the  information  needed  for 
the  accurate  estimation  of  the  fluxes  is  local,  and  that  the  use  of  local  estimates  does  not 
lead  to  variations  which  are  much  greater  than  those  observed  with  global  estimates.  This 
is  consistent  with  the  observations  of  Langford  [30]  in  spectral  based  optimal  LES. 

With  these  observations  about  local  estimates  in  mind,  one  can  make  suggestions  about 
the  stencils  to  be  used  in  actual  simulations.  For  the  viscous  fluxes  a  standard  linear  model  is 
adequate,  while  a  quadratic  model  (as  expected)  must  be  used  for  the  convective  fluxes.  To 
improve  the  accuracy  of  the  models,  the  extension  of  the  stencil  should  be  primarily  in  the 
direction  normal  to  the  face.  If  further  accuracy  is  desired,  the  stencil  can  be  widened  and, 
in  this  case  only,  nonstandard  models  should  be  used.  This  last  recommendation  appears 
more  important  for  collocated  grids  than  for  staggered  ones. 

C.  A  Posteriori  Analysis  of  Local  Estimates 

The  choice  of  local  models  for  the  a  posteriori  analysis  is  based  on  the  results  from  the 
previous  section.  The  models  analyzed  a  posteriori  are  listed  in  Table  II  with  their  a-priori 
variations  displayed  in  Table  III.  The  nomenclature  of  the  models  is  the  following:  the  first 
letter  indicates  whether  the  model  is  for  a  staggered  or  collocated  grid,  the  second  character 
(2,  4,  6  or  w)  indicates  the  stencil  geometry  and  the  third  letter,  if  present,  indicates  the 
use  of  standard  (s)  or  nonstandard  (n)  finite-volume  events  in  the  model.  Model  descriptors 
without  a  third- character  refer  to  models  using  only  standard  events.  Note  that  the  stencils 
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are  only  given  for  four  fluxes  shown  in  Table  II  because  by  isotropy,  these  four  fluxes  can  be 
rotated  to  define  the  stencil  for  any  of  the  required  fluxes  (see  Section  II  A).  For  this  reason, 
there  is  no  v-dependency  in  any  of  the  standard  models  shown. 

In  a  priori  analysis,  lengthening  the  model’s  stencil  yielded  the  largest  improvements 
in  the  accuracy  of  the  flux.  Widening  the  stencil  did  not  bring  such  benefits,  unless  non¬ 
standard  terms  were  used.  With  these  characteristics  in  mind,  most  of  the  models  for  a 
posteriori  analysis  were  constructed  with  long,  thin  stencils.  Wider  models  (models  Cws, 
Cwn,  Sws  and  Swn)  were  used  to  verify  if  extending  the  stencil  in  the  direction  parallel  to 
the  face  was  indeed  ineffective  for  error  reduction.  Most  of  the  a  posteriori  models  only  con¬ 
tain  standard  terms,  but  four  pairs  of  models  (model  pairs  C6s-C6n,  Cws-Cwn,  S6s-S6n  and 
Sws-Swn),  with  identical  stencil  geometries,  were  simulated  with  standard-only  and  stan- 
dard+nonstandard  terms.  These  results  allow  the  investigation  of  the  effect  of  nonstandard 
terms  on  model  accuracy. 

The  a  posteriori  results  were  obtained  using  a  third-order  Runge-Kutta  scheme  for  time 
integration.  To  maintain  turbulence,  a  negative  viscosity  forcing  term  was  used  to  stir  the 
flow.  The  forcing  was  identical  to  the  one  used  in  the  original  DNS  [10].  The  LES  statistics 
were  collected  for  approximately  60i,  with  samplings  separated  by  0.15£,  where  t  is  the 
eddy-turnover  time  defined  by 

t-£  (39) 

2e 

The  initial  fields  for  the  simulation  were  obtained  by  filtering  DNS  fields. 

Pressure  was  determined  by  projecting  the  velocity  field  onto  a  discretely  divergence-free 
space,  using  an  approximate  divergence  operator.  As  pointed  out  by  Langford  &  Moser 
[27],  such  a  projection  introduces  a  modeling  error  because  the  filtered  field  does  not  exactly 
satisfy  any  continuity  constraint.  As  suggested  by  their  results,  a  two  or  three-point  stencil 
(for  staggered  and  collocated  grids,  respectively)  for  the  approximate  divergence  was  used, 
which  provides  an  accurate  approximation,  at  least  in  the  staggered  grid  case.  By  symme¬ 
try,  these  approximate  divergence  operators  are  equivalent  to  the  second-order  divergence 
approximations  commonly  used  in  finite-volume  methods.  Unfortunately,  many  of  the  col¬ 
located  models  were  not  numerically  stable  and  results  were  obtained  only  for  few  of  them 
(models  C2,  C4  and  Cwn).  The  reason  for  this  instability  appears  to  be  connected  to  the 
treatment  of  pressure,  leading  to  the  commonly  observed  odd-even  decoupling  in  collocated 
models  without  enough  dissipation  at  the  highest  wavenumbers. 
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The  three-dimensional  energy  spectrum  is  shown  in  Fig.  3.  Most  of  the  models  agree  very 
well  with  the  filtered  DNS  energy  spectrum,  with  small  errors  in  the  high  wavenumber  region. 
The  smallest  stencils,  represented  by  models  C2  and  S2,  do  not  seem  to  incorporate  enough 
information  to  represent  the  higher  wavenumbers.  The  collocated  model  C2  is  especially 
inaccurate,  only  wavenumbers  up  to  around  7  are  correctly  represented.  The  staggered 
model  S2  is  better,  but  there  seems  to  be  excess  dissipation  at  the  higher  wavenumbers. 
The  larger  staggered  stencils  are  remarkably  better.  With  model  S4,  the  LES  spectrum  is 
accurate  throughout  the  wavenumber  range,  and  model  S6  improves  upon  those  results  only 
slightly.  The  collocated  C4  model  only  improves  slightly  upon  the  C2  model,  and  is  not 
nearly  as  accurate  as  the  staggered  models. 

Surprisingly,  widening  the  stencil,  as  in  model  Sws,  caused  a  large  discrepancy  in  the 
energy  of  the  higher  wavenumbers.  Inspection  of  the  stencils  for  the  Sws  model  showed  that 
the  linear  terms  for  the  Muw  fluxes  had  an  anti-diffusive  character  for  large  wavenumbers. 
This  is  likely  to  lead  to  poor  large-wavenumber  performance,  and  may  lead  to  numerical 
instabilities.  There  is  no  a-priori  guarantee  that  the  optimal  estimation  procedure  will 
produce  numerical  schemes  that  are  stable.  Such  stability  constraints  could  be  included  in 
the  optimization  and  this  should  be  explored  in  future  research. 

The  models  with  nonstandard  terms  (Cwn,  S6n  and  Swn)  were  markedly  inaccurate  when 
compared  with  the  similar  models  which  used  only  standard  terms  (S6s  and  Sws),  especially 
at  high  wavenumbers.  The  value  of  nonstandard  terms  in  reducing  the  a  priori  error  is 
dependent  on  the  correlations  between  different  velocity  components.  Apparently  in  the 
LES,  these  cross  correlations  were  not  maintained  causing  the  nonstandard  terms  to  degrade, 
rather  than  enhance  performance.  An  examination  of  the  coefficients  for  the  nonstandard 
terms  reveals  that  many  of  them  involve  what  appear  to  be  discrete  approximations  to 
du/dx,  dv/dy  and/or  dw/dz,  so  as  to  suggest  that  the  quantities  being  modeled  can  be 
rewritten  using  continuity.  Thus,  it  may  be  that  modeling  errors  inherent  to  the  continuity 
representation  contribute  to  the  poor  performance  of  nonstandard  models. 

The  energy  transfer  by  the  models  can  be  analyzed  by  computing  the  third-order  structure 
function 

S3(r)  =  ((u(x  +  r)  -u(x)f),  (40) 

which  is  shown  in  Fig.  4.  The  structure  function  results  for  the  smaller  models  (C2  and 
S2)  are  not  in  agreement  with  the  DNS  data  for  small  separations.  This  reveals  that  the 
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energy  transfer  is  not  adequately  represented  by  stencils  smaller  than  4  in  the  direction 
normal  to  the  face,  and  explains  the  deficiencies  encountered  in  the  energy  spectrum  of  the 
smaller  models.  The  remaining  staggered  models  reproduce  S3  quite  accurately,  while  the 
C4  collocated  model  underpredicts  it  at  all  separations  and  the  Cwn  collocated  model  is 
quite  inaccurate  for  small  separations.  For  the  staggered  grids,  the  widest  stencil  (model 
Sws)  is  slightly  better  than  the  thin  stencils  (models  S4  and  S4s),  but  the  difference  is  almost 
negligible.  The  nonstandard  models  are  again  less  accurate  than  their  standard  counterparts. 

In  Appendix  B,  the  model  coefficients  for  the  S2  and  S4  models  are  given.  The  S4  model 
is  given  because  it  is  the  smallest  stencil  that  produces  good  agreement  with  the  filtered 
DNS  data,  for  both  the  spectra  and  the  structure  function.  The  S2  model  is  shown  for 
comparison  with  a  standard  second-order  finite-volume  scheme.  One  important  difference 
between  the  S2  model  and  the  standard  second-order  staggered  approximation  is  that  the 
optimal  model  is  more  dissipative,  primarily  due  to  the  linear  term  in  the  convective  fluxes. 
This  is  expected  because  the  optimal  model  includes  the  effects  of  the  subgrid  turbulence, 
which  is  primarily  dissipative.  Also  shown  in  Appendix  B  is  the  procedure  for  constructing 
the  S2  model,  that  can  be  easily  extended  to  all  other  optimal  models  discussed  in  this 
paper. 

In  Figs.  5  and  6,  results  for  the  optimal  finite-volume  LES  models  are  compared  with  re¬ 
sults  obtained  with  the  dynamic  Smagorinsky  subgrid  scale  model.  The  dynamic  Smagorin- 
sky  results  were  obtained  on  staggered  grids  using  the  same  code  used  for  the  optimal  finite- 
volume  LES  cases,  employing  second-  and  fourth-order  accurate  discretization  schemes  for 
the  fluxes  and  subgrid  models.  In  Fig.  5,  the  dynamic  Smagorinsky  model  dampens  the 
high- wavenumber  spectra  excessively  when  compared  to  the  optimal  models  and  the  filtered 
DNS  results,  while  the  spectra  for  optimal  and  dynamic  Smagorinsky  models  is  comparable 
up  to  wavenumbers  around  18. 

In  Fig.  6,  the  structure  function  results  indicate  that  the  smaller  stencils,  S2  and  the 
second-order  accurate  dynamic  Smagorinsky,  are  not  appropriate  for  an  accurate  represen¬ 
tation  of  the  energy  transfer  in  this  LES.  For  these  small  stencils,  the  dynamic  Smagorinsky 
model  is  better  than  the  optimal  model  at  small  separations,  while  overpredicting  the  struc¬ 
ture  function  by  a  large  margin  when  the  separation  is  larger.  When  using  the  larger  stencils, 
as  in  the  S4  optimal  model  and  the  fourth-order  accurate  dynamic  Smagorinsky  model,  the 
energy  transfer  of  the  fourth-order  dynamic  Smagorinsky  model  is  slightly  more  accurate 
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than  that  of  the  second-order  model,  but  the  structure  function  is  still  overpredicted  in  most 
of  the  separation  range.  On  the  other  hand,  the  optimal  model  S4  structure  function  is  very 
close  to  the  filtered  DNS  results  for  all  separations.  Therefore,  for  this  example,  the  optimal 
model  S4  is  more  accurate  than  the  dynamic  Smagorinsky  subgrid  model  using  a  standard 
finite-volume  discretization  scheme  with  an  identically  sized  stencil. 

IV.  DISCUSSION 

An  optimal  coarse  finite-volume  LES  modeling  technique  has  been  developed  and  ap¬ 
plied  to  isotropic  turbulence.  This  approach  addresses  some  of  the  issues  encountered  in 
current  LES  modeling,  namely:  ambiguous  definition  of  filtering  and  treatment  of  numer¬ 
ical  discretization  errors.  While  many  current  LES  models  allow  the  filter  to  be  defined 
implicitly  by  the  numerics,  the  optimal  finite-volume  LES  technique  rigorously  defines  the 
filter  (equation  (9))  and  the  resulting  model  is  strongly  dependent  on  this  choice,  as  it  must 
be.  Numerical  discretization  issues  usually  encountered  in  coarse  representations  of  LES 
are  resolved  in  the  optimal  technique  by  handling  the  Navier-Stokes,  discretization  and  sub¬ 
grid  effects  simultaneously.  In  this  way,  the  combined  numerical  and  modeling  errors  are 
minimized. 

The  results  obtained  while  testing  the  optimal  finite-volume  staggered  models  are  encour¬ 
aging.  The  energy  spectra  and  third-order  structure  functions  computed  a  posteriori  are  of 
comparable  accuracy  to  those  obtained  by  Langford  [30]  using  a  spectral  numerical  method, 
a  sharp  Fourier-cutoff  filter  and  an  optimal  model  for  the  subgrid  stress.  The  spectral  nu¬ 
merical  method  used  in  Langford  [30]  avoids  the  discretization  errors  usually  encountered  in 
LES.  The  current  approach  achieves  similar  results,  indicating  that  the  numerical  limitations 
one  expects  in  a  finite-volume  representation  have  largely  been  avoided.  The  a-posteriori 
results  for  the  optimal  model  were  more  accurate  than  those  obtained  with  the  dynamic 
Smagorinsky  model  using  either  standard  second-  or  fourth-order  discretizations  for  the  flux 
terms.  The  dynamic  Smagorinsky  results  show  high-wavenumber  discrepancies  in  the  en¬ 
ergy  spectra  that  are  not  observed  in  the  data  from  the  optimal  models.  This  behavior  is 
consistent  with  the  inability  of  the  coarse  low-order  standard  numerical  method  to  resolve 
the  high-wavenumber  components  of  the  resolved  field. 

While  the  performance  of  the  optimal  finite  volume  LES  evaluated  here  was  quite  promis- 
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ing,  there  are  a  few  issues  that  need  to  be  addressed.  In  a  few  isolated  cases,  the  optimal 
estimation  procedure  yielded  model  coefficients  with  poor  numerical  properties,  and  in  some 
cases  the  models  were  unstable.  Methods  that  constrain  the  optimal  procedure  to  produce 
only  stable  numerical  schemes  should  be  pursued  in  future  research.  Instabilities  also  arose 
in  collocated  models,  probably  due  to  the  implementation  of  the  continuity  constraint.  Tra¬ 
ditional  finite- volume  schemes  have  similar  difficulties  that  commonly  lead  to  checkerboard 
decoupling  of  the  pressure.  In  future  research  with  collocated  grids,  stabilization  procedures 
for  computation  of  the  pressure  (e.g.  Rhie  &  Chow  [36])  may  be  useful. 

The  modeling  procedure  described  here,  in  which  DNS  statistical  data  is  used  as  input 
to  the  optimal  modeling  formulation,  is  clearly  not  useful  for  generating  practical  LES  mod¬ 
els,  unless  the  need  for  DNS  data  can  be  eliminated.  Because  the  results  of  the  current 
study  indicates  that  our  approach  can  produce  accurate  LES,  such  a  generalization  is  being 
pursued.  A  promising  approach  to  eliminating  the  need  for  LES  data  is  to  evaluate  their 
required  statistical  quantities  theoretically.  The  correlations  required  as  input  to  the  mod¬ 
eling  (see  Appendix  B)  are  surface  and  volume  integrals  of  multi-point  second,  third  and 
fourth  order  velocity  correlations.  Using  a  combination  of  Kolmogorov  theory  (the  |  and 
|  laws  for  the  second  and  third  order  correlations  respectively),  the  quasi-normal  approx- 
imation  (fourth  order  correlations)  and  a  generalized  dynamic  procedure,  it  appears  that 
the  required  correlations  can  be  determined  under  the  assumption  that  filter  width  is  in 
an  inertial  range,  and  that  the  small  scales  are  isotropic.  Once  these  multi-point  velocity 
correlations  are  determined,  they  can  be  integrated  as  needed  for  the  finite  volume  grid 
being  used  in  any  particular  problem,  and  the  optimal  modeling  coefficients  specific  to  that 
grid  can  then  be  found.  We  envision  this  as  a  preprocessing  step  for  a  simulation,  similar 
to  constructing  and  storing  the  stencil  coefficients  for  a  standard  differencing  scheme  on  a 
general  mesh.  When  the  assumptions  of  an  inertial  range  and  small  scale  isotropy  are  not 
valid  (e.g.  near  walls),  then  further  modeling  will  be  required.  However,  it  is  still  a  matter 
of  modeling  the  small-separation  multi-point  velocity  correlations. 

It  is  useful  to  contrast  the  development  of  the  finite  volume  optimal  LES  formulation 
described  here  with  common  finite  volume  numerical  schemes  and  to  the  MILES  approach 
introduced  by  Boris,  et  al.  [37].  In  the  current  formulation,  the  representation  in  terms  of 
volume  averaged  velocities  is  identical  to  standard  finite  volume  techniques.  But,  because 
the  volume  size  is  by  definition  large  compared  to  the  smallest  scale  of  motion,  it  is  not  valid 
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to  approximate  fluxes  using  standard  reconstruction  techniques,  which  are  based  on  the 
assumption  that  the  velocity  is  smooth  on  the  scale  of  the  volume  size.  Another  circumstance 
in  which  this  assumption  is  violated  is  shock  capturing.  In  this  instance,  special  numerical 
approximations  have  been  devised  to  account  for  the  sub-grid  features  of  the  flow  (the 
shocks),  based  on  knowledge  of  the  properties  of  shocks.  Thus  shock  capturing  schemes 
are  designed  to  recover  the  shock  jump  conditions,  keep  the  numerical  shock  as  thin  as 
possible  on  the  grid  (of  order  the  grid  spacing)  and  preserve  the  monotonicity  of  the  shock 
solution.  Subgrid  turbulence  is  more  complex  and  our  knowledge  of  the  properties  of  subgrid 
turbulence  is  only  statistical.  It  is  thus  this  statistical  information  that  is  used  in  the  optimal 
LES  formulation.  The  MILES  formulation  first  introduced  by  Boris,  et  al.  [37]  is  based  on 
the  observation  that  one  of  the  statistical  properties  of  subgrid  turbulence  is  that,  like  shocks, 
it  is  dissipative.  Numerical  schemes  similar  to  those  designed  for  shock  capturing  are  thus 
used  in  MILES,  with  the  dissipative  properties  adjusted  to  more  closely  resemble  those  of 
turbulence.  However,  in  general,  subgrid  turbulence  has  numerous  statistical  properties 
that  impact  large-scale  dynamics,  besides  dissipation.  In  particular,  with  finite  volume 
representations  like  those  described  here,  the  subgrid  affects  every  aspect  of  the  filtered 
dynamics.  The  optimal  approach  allows  models  to  be  constructed  that  reproduce  most  of 
the  dynamically  significant  statistical  properties  o  priori  [11]. 

The  finite  volume  optimal  LES  formulation  described  here  was  pursed  rather  than  the 
theoretically  more  straight-forward  formulations  based  on  spectral  representations  because 
a  spectral  formulation  cannot  easily  be  applied  in  flows  with  complex  geometries.  While,  the 
current  study  involves  only  a  very  simple  flow  (isotropic  turbulence)  on  a  uniform  Carte¬ 
sian  grid,  the  approach  is  clearly  applicable  in  general  geometries  on  general  grids.  The 
only  limitation  is  the  availability  of  the  required  statistical  data.  In  this  study,  statistical 
data  was  obtained  from  DNS,  but  current  research  is  directed  at  providing  this  information 
theoretically. 
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APPENDIX  A:  ANALYSIS  OF  ESTIMATES  BASED  ON  SINGLE  AND  DUAL¬ 
FLUX  VARIATIONS 

While  the  dual-flux  variation  d\\  is  more  relevant  to  the  LES  evolution,  using  this  variation 
in  the  estimation  procedure  is  not  without  problems,  as  demonstrated  by  the  following  one¬ 
dimensional  model  problem.  Notation  for  this  problem  is  illustrated  in  Fig.  7 .  The  variable 
w  is  a  scalar  and  M  is  the  x-component  of  a  vector  field. 

The  model  system  is  considered  to  be  homogeneous,  with  a  reflectional  symmetry  in  its 
single  spatial  direction.  It  follows  that  the  symmetries  of  this  system  are 


(' w+w+ )  =  {ww) , 

(Al) 

(w+w)  =  (ww~) , 

(A2) 

-  (M+w)  =  (M+w+) . 

(A3) 

The  exact  evolution  of  the  model  system  is  w  =  M+  —  M~ ,  where  w  = 

dw/dt.  Therefore, 

the  mean  dissipation  e  in  this  system  is 

e  =  {ww) 

=  (M+w)-(M~w) 

=  2  (M+w) 

=  -2(M+w+). 

(A4) 

1.  Estimate  based  on  Single-Flux  Variation  (d|) 

If  the  linear  stochastic  model  of  M+  given  event  data  w  and  w+  is: 

m+  =  aw+  +  (3w . 

(A5) 

The  coefficients  a  and  /?  can  be  determined  by  the  stochastic  estimation  procedure.  When 
minimizing  the  single-flux  variation  d\  =  M+  —  m+,  the  stochastic  estimation  procedure 
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yields  the  following  linear  system 


(io+w+)  (ww+)  \  I  a 
(w+w)  (ww)  I  \  (3 


(. M+w+ ) 

( M+w ) 


(A6) 


FYom  the  symmetries  (A1-A3)  it  follows  that  the  solution  for  the  coefficients  is 

(M+w+) 

a  =  —  p  = 


(A7) 


(ww)  —  (ww+) 

The  estimated  evolution  of  the  system  is  w  —  m+  —  m~ ,  so  that  the  estimated  dissipation 


e  is 


e  =  (^ww^ 


—  ( m+w )  —  ( m~w ) 

=  a  (( w+w )  —  (ww))  —  a  ((ww)  —  ( w~w )) 
=  2a  ((ww+)  —  (ww))  . 


(A8) 


Substituting  equation  (A7)  into  equation  (A8)  and  comparing  with  the  expression  given  by 
(A4),  is  is  clear  that  e  =  e. 

The  shortcoming  to  this  method  of  estimation  is  the  possibility  of  fluctuations  of  M+ 
that  exactly  cancel  fluctuations  of  M~.  Unfortunately,  part  of  the  “effort”  of  computing 
the  estimate  goes  to  estimating  the  part  that  cancels.  An  improvement  may  be  achieved  by 
minimizing  the  variation  in  the  dual-flux,  M+  -  M~,  rather  than  the  single-flux  variation 
in  only  M+. 


2.  Estimates  based  on  the  Dual-Flux  variation  (dy) 

When  using  the  dual-flux  variation,  the  estimate  is  of  the  following  form 

m+  —  m~  =  aw+  +  (3w  —  aw  —  /3w~  (A9) 

=  aw+  +  (/?  —  a)w  —  pw~  (A10) 

=  a(w+  —  w)  - f  P(w  —  w~).  (All) 

This  can  be  regarded  as  a  linear  stochastic  model  of  (m+  —  m~)  given  event  data  (w+  —  w) 
and  (w  —  w~). 
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The  estimation  equations  in  this  case  are  given  by  the  linear  system 


{{w+  —  w)(w+  —  w))  ((w  —  w  ){w+ —  w))  \  (  a 
{{w+  —  w)(w  —  w~))  {(w  —  w~)(w  —  w~))  J 

((M+  —  M~)(w+  —  w)} 
(( M+  —  M~){w  —  w~)) 


The  solution  for  this  system  is 
a  =  ~P  — 


{{M+  -  M~)(w+-w)) 


(( w  —  w~)(w  —  w~ ))  —  (( w  —  w  )(w+  —  w ))  ’ 
and  the  estimated  dissipation  is 

e  =  (ww^ 

=  ((m+  —  mT)w) 

=  a  (((w+  —  w)w)  —  (( w  —  w~)w )) 

=  a  ((ww+)  —  2  (ww)  +  (ww~)) 

=  2a  (( wu>+ )  —  {ww})  . 


(A12) 


(A13) 


(A14) 


Note  that  (A14)  is  identical  to  equation  (A8)  of  the  single-flux  estimate.  The  value  of  a  for 
the  previous  solution,  equation  (A7),  was  such  that  e  =  e.  The  value  of  a  for  the  current 
technique,  equation  (A13),  is  different.  The  current  solution  for  a  given  in  (A13)  depends  on 
(w~w+),  while  the  previous  solution,  equation  (A7),  does  not.  It  follows  that  the  coefficients 
a  for  the  two  cases  are  in  general  different,  and  that  for  this  method,  d  ^  d. 

Therefore,  a  dual-flux  estimate  written  in  conservative  form  will  not  correctly  match 
dissipation.  One  can,  however,  choose  to  estimate  the  dual-flux  in  a  nonconservative  form, 
i.e. 


m+  —  m  =  aw+  +  pw  +  'yw  . 


In  this  case,  the  estimation  equations  are 


^  (w+t u~)  { ww~ )  { w~w~ )  ^ 
(w+w)  {ww)  {w~w) 
y  {w+w+}  {ww+)  { w~w+ )  j 


P 

w 


(  {{M+  -  M~)w+)  } 
{{M+  -  M~)w) 
y  ((M+  -  M~)w~)  j 


(A15) 


(A16) 
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Symmetries  dictate  that  a  =  7,  but  nothing  is  implied  about  the  relationship  to  (3.  Without 
actually  solving  these  equations,  we  will  simply  consider  the  second  row  in  (A16).  It  follows 
that 

((ra+ —  m~)w)  =  ((M+  —  M~)w) 

d  =  d.  (A17) 

This  estimate  does  not  suffer  from  cancellation  problems,  and  it  has  the  correct  dissipation. 
However,  it  is  not  conservative. 

In  the  three  cases  shown  above,  it  was  demonstrated  that  there  are  several  choices  for 
local  estimation  techniques,  but  they  all  suffer  from  a  deficiency.  One  can  construct  a  con¬ 
servative  single-flux  estimate  that  correctly  matches  dissipation,  but  some  of  the  estimation 
power  goes  to  minimize  an  error  in  portions  of  the  flux  that  trivially  cancel.  Or  one  can 
construct  a  conservative  dual-flux  estimate  that  is  aware  of  the  cancellations,  but  dissipa¬ 
tion  is  not  correctly  matched.  Finally,  one  can  construct  a  dual- flux  estimate  that  minimizes 
the  most  relevant  error  and  correctly  matches  dissipation,  but  the  resulting  estimate  is  not 
conservative. 

There  is  yet  another  alternative.  One  can  construct  a  conservative  dual-flux  estimate, 
subject  to  the  constraint  that  dissipation  is  matched,  i.e.  that 

a  ((u>+  —  w)w )  +  (3  {{w  —  w~)w )  =  ((M+  —  M~)w )  .  (A18) 

For  the  simple  example  explored  here,  the  procedure  reduces  to  the  conservative  single-flux 
estimate.  However,  as  additional  nonlocal  velocity  data  is  incorporated  into  the  estimate, 
the  procedure  yields  a  distinct,  fourth  estimation  procedure.  Also  note  that  in  the  limit 
of  a  global  estimate,  all  four  estimation  procedures  are  identical.  This  can  be  understood 
by  noting  that  for  the  global  estimates,  one  can  regard  the  procedure  as  the  simultaneous 
estimation  of  fluxes  at  all  locations,  not  just  flux  at  a  single  location.  Then,  the  combination 
of  fluxes  at  opposite  faces  is  merely  a  linear  transformation  of  the  quantity  being  estimated. 


APPENDIX  B:  MODELS  S2  AND  S4  FOR  STAGGERED  GRIDS 

The  coefficients  for  the  optimal  finite-volume  LES  models  S2  and  S4  as  formulated  in 
the  equations  below,  are  presented  in  table  IV,  with  coefficients  for  the  standard  second- 
order  finite-volume  scheme  shown  for  comparison.  The  optimal  models  are  strictly  valid 
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only  for  the  test  case  presented  (microscale  Reynolds  number  Re  a  =  164,  Taylor  microscale 
A  =  0.209,  kinetic  energy  <j2/ 2  =  41.1,  dissipation  e  =  62.9  and  grid  spacing  A  =  7r/16). 

+  C(Uij_ltk  +  +  u  +  ut,j+i,k)  +  D(  'U'iJ- 2,fc  "t”  ^i,j,k- 2  ”1“  ”1“  ^i,i+2,fc) 

+  fc)  +  -  wi+1>Mwi+ltiifc)  +  G(u^lijjkuitj>k  -  uiJikui+1Jtk ) 

+  /U'itj—l,k'Vitjtk  "I”  l^i  —  l,jtk  ^i,j,k—l^ifj,k 

l,j+l,fc 

+  +  VijtkUi)j)k  — 

+  Wi-xtjikUijyk  —  +  Wi,j,k'U'itjyk  ~~  ^i^.fc+i^i^fc+l)  ”  Pa; 

dv  A/  x  /  . 

_  =  )i-1>fc  +  ^>i+1,fc)  +  5(^-2,*  +  ^i,i+2,fc) 

+  C(Vi_ltjfk  +  Vitjlk- 1  +  +  ^i+l,j,fc)  +  -0(^i^2,i,fc  +  Vi,j,k- 2  +  ^i,j,fc+2  +  ^i+2,i,fc) 

“h  E(Vi +  F{yiyj-\ykViyj„xyk  Vi,j+l,kVi,j+l,k)  “1“  G(Vij-iykVitjyk  Vitj,kVi,j+ltk) 

+  H (—Vij-Z'kVij-i^  +  vitj+itkVij+2tk)  H"  I(vi,j-2  2,fc  ^i,jf+2,  fc^i,j+2,fc) 

H”  ^iV'iyj-ltk'^iJtk  “1”  'U'i,jikViijik  'Ui+i)j-ltk'Vi+l}jyk  'U'i+l,j,k'Ui+ltjtk 

“f“  “1“  “I”  ^itjyk- 1  ,fc  ”1”  ^i,j,k 
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~fa  =  1  +  W«,fc+l)  +  £(w<J,fc-2  +  W^>fc+2) 

+  G{wi_ltjik  “{“  Wij—  iffc  “f*  ^i,jr+l,fc  +  ^<+ i,j,k)  +  D(w{_  2tj>k  +  witj_2>k  +  Wi(J+2,fc  +  Wi+^j.k) 

“H  E{y^i,j,k)  -f"  F(Wijik-iWitj'k-i  G{^wijk_iWijk  WijikWijik+i) 

+  +  wiijM1  witjik+2)  +  I(wiJik_2witj'k_2  -  wiijtk+2wiJik+2) 

J (Ui tjtk  —  l'Witjtk  “I”  'U'i+l,jtk—  l'Mi+l,j,k  'U'i+l1jtk'Wi+ltj,k 

+  Vij'k-rWij'k  +  viJikwiijtk  ~  Viij+ltk.1Wi>j+lik  -  Vitj+1  ,kwiJ+hk 
^1—1, "t”  l,j,k'U’i,j,k  “b  l,fc^s,j,fc—  1  4"* 

'Mitj,k'U/i  +  l,j,k—l  ^  itj  ,k^i+\ ,j,k  lMi)jtk'Uilj+ltk  —  l  tUiJ  1  ,fc  )  Pz 


The  procedure  for  computing  model  S2  is  described  below,  the  extension  to  the  other 
models  is  straight  forward.  Each  of  the  basic  fluxes  (Muw,  Mww,  Mgu/dz  and  Mgw/gz)  must 
obey  the  stochastic  estimation  equations  (23)  and  (24),  which  are  repeated  here  for  conve¬ 
nience: 


(AT)  =  (mi>  (Bl) 

<(Mi  -  m^Ef)  =  0  (B2) 

Note  that  the  mean-preserving  condition  (Bl)  may  be  incorporated  into  (B2)  by  adding 
a  constant  term  to  the  event  vector  Ex  and,  therefore,  the  estimation  equations  can  be 
expressed  as: 

(mEt)  =  (MEt)  (B3) 

The  i  superscripts  may  be  removed  since  the  correlations  are  averaged  over  all  volumes.  For 
clarity’s  sake,  the  following  notation  will  be  used  for  the  filtered  velocity  components: 


U-  —  u  (x) 

u+  =  u(x  - f  (0, 0,  A)) 

W-  =  w  ( x ) 

w+  =  w  (x  +  (0, 0,  A)) 


(B4) 

(B5) 

(B6) 

(B7) 

(BS) 

(B9) 
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The  event  vectors  for  each  basic  flux  are: 

M UW 


Mdu/dz : 


Mdw/dz  • 


E  = 


E  = 


(  i  \ 

U-Wr 

U-VJi 

u+wr 

\U+Wi  f 

(  1  \ 

W- 

w+ 

VJ-W- 

w+w+ 
\U!-W+  J 


E  = 


E 


U- 

w 


W 


(BIO) 


(Bll) 


(B12) 


(B13) 


Note  that  the  equation  for  the  model  (18)  may  be  written  in  the  following  form  (summa¬ 
tion  on  repeated  indices): 

m*  =  CkEl  (B14) 

where  C  is  a  vector  of  coefficients.  Therefore,  the  estimation  equation  (B3)  may  be  expressed 
as: 

Ck  (EkEt)  =  (MEi)  (B15) 
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or  in  matrix  notation: 


(ExEi)  (E1E2)  . 
(EtEi)  (E2E2)  . 

•  •  (EiEn)~ 
..  {E2En) 

'Ci 

C2 

— 

‘{MEi) 

< ME2 ) 

(EnEi)  (EnE2)  . 

•  •  (EnEn)_ 

Pn. 

_(MEn)_ 

(B16) 


where  n  is  the  number  of  events  for  the  model.  Notice  that  the  correlation  matrix  is  sym¬ 


metric. 

Therefore,  for  model  S2,  the  matrix  form  of  the  estimation  equations  is  (since  the  matrices 
are  symmetric,  only  the  upper  triangle  of  each  is  shown): 

Muw: 


1  (u_)  («+)  (u-w r)  ( U-Wl )  (u+W  r)  { U+Wt ) 

qMuw 

{Muw) 

(u2_)  ( U-U+ )  ( U2_Wr )  (u2_Wi)  (U-U+Wr)  { U-U+Wi ) 

{MUWU— ) 

(u2+)  (u+u-wr)  (u+u-.wi)  (u\wr)  (■ u\wi ) 

qMu-w 

{Muwu+  ^ 

{(U-Wrf)  (u2_WrWi)  (U-U+W2)  (U-WrU+Wi) 

qMuw 

= 

(Muwu-Wr) 

(( U-Wi )2)  ( U-WiU+Wr )  {U-U+wf) 

qMuw 

(MuwU-Wi ) 

((u+wr)2)  ( u\wrwi ) 

qMuw 

(. Muwu+Wr ) 

((u+wi)2)  _ 

_(MUWU+Wi)_ 

1  (w-)  (w+)  (w2_)  (wl)  {' W-W+ ) 

M u>  u? 

(Mww) 

(w2_)  ( W-W+ )  (wt)  {w-w\)  (w2_w+) 

qMww 

(. Mwww- ) 

(w2+)  (w+w2)  (wl)  (w-wl) 

(Mwww+) 

(wt)  ((w-w+)2)  (wtw+) 

(Mwww2_) 

(wl)  (w-wl) 

A^ww 

(Mwww+) 

((w.w+)2)_ 

(MWVJw-w+)_ 

Mdu/dz  • 


1  (u_)  (u+) 

qMquiqz 

(Mdi ijdz) 

(u2_)  (■ u+u _) 

^i^du/dz 

= 

( Mdu/dzU -) 

to) . 

/~iMgufQz 

°3 

(Mdu/dz'Ui+) 

(B18) 


(B19) 
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Mdw/dz  * 

1  (u>_)  (w+)  1  C™dw/dz  {Mdw/dz) 

{w2_)  { w+w _}  Q^w/az  _  {Mdw/dzwJ)  (B20) 

K)  cfW_  (Mdw/dzw+)_ 

The  solution  for  these  estimation  systems  gives  us  the  models  for  the  basic  fluxes  corre¬ 
sponding  to  model  S2.  Many  of  the  correlations  appearing  in  the  estimation  equations  may 
be  simplified  by  using  homogeneity  and  isotropy,  such  as: 

(«+)  =  («->  =  (™+)  =  ( V -)  =  0 

«>  =  <-2->  =  «>  -  <«£> 

( u+ut )  =  -  («-«+) 

The  extension  of  the  estimation  procedure  to  other  model  inputs  and  stencils  is  trivial, 
since  the  estimation  equations  Eire  only  based  on  the  event  vector  E.  When  constructing 
estimation  equations  for  a  different  model,  all  one  must  do  is  form  an  event  vector  with  all 
model  inputs,  compute  the  correlation  matrix  ( E{Ej )  and  the  right-hand  side  vector  { MEi ), 
and  solve  the  resulting  linear  system  (( EiEj )  Cj  =  {MEi))  for  the  coefficients  C.  In  this 
work,  Gaussian  elimination  was  used  to  solve  the  linear  system. 
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TABLE  I:  Relative  estimation  variations  for  various  finite-volume  fluxes.  Variations  are  shown  for 
linear  estimates,  for  estimates  with  local  quadratic  products  (LQ),  and  for  estimates  with  nonlocal 
quadratic  products  (NLQ)  formed  from  velocities  in  cells  separated  by  a  distance  of  up  to  one  cell 
width.  Also  shown  are  simplified  estimates  (f),  for  which  only  terms  from  standard  finite- volume 


approximations  are  included. 

grid  arrangement 

relative 

variation 

and  estimation  type 

M-WVJ 

Mdu/dz 

Mdw/dz 

collocated  linear 

0.9983 

0.9983 

0.4725 

0.3170 

collocated  linear  f 

0.5065 

0.3554 

collocated  LQ 

0.4943 

0.4530 

0.4720 

0.3165 

collocated  LQ  f 

0.5341 

0.4716 

collocated  NLQ 

0.4430 

0.3530 

0.4660 

0.3137 

collocated  NLQ  f 

0.5040 

0.4447 

staggered  linear 

0.9983 

0.9983 

0.4622 

0.3296 

staggered  linear  f 

0.5065 

0.3554 

staggered  NLQ 

0.4802 

0.4174 

0.4588 

0.3275 

staggered  NLQ  f 

0.5033 

0.4494 

TABLE  II:  Stencils  for  local  models  used  in  the  a  posteriori  analysis,  f  Models  with  standard 
terms  only.  *  Models  which  diverged  during  simulation. 


model 

grid  type 

flux 

u 

Stencil  dimensions 

V 

w 

C2  f 

collocated 

all 

1x1x2 

- 

1x1x2 

C4  f 

collocated 

all 

1x1x4 

- 

1x1x4 

C6s  f  * 

collocated 

all 

1x1x6 

- 

1x1x6 

Cws  f  * 

collocated 

all 

3x3x4 

- 

3x3x4 

C6n  * 

collocated 

all 

1x1x6 

1x1x6 

1x1x6 

Cwn 

collocated 

all 

3x3x4 

3x3x4 

3x3x4 

MUW)  Mdu/dz 

1x1x2 

- 

2x1x1 

S2  f 

staggered 

MWW 5  Mdw/dz 

- 

- 

1x1x2 

MUW)  Mdu/dz 

1x1x4 

- 

2x1x1 

S4f 

staggered 

Mwwj  Mdw/dz 

- 

- 

1x1x4 

J  M-du/dz 

1x1x6 

- 

2x1x1 

S6sf 

staggered 

)  M. dw/dz 

- 

- 

1x1x6 

Muw ,  Mdu  /  dz 

3x3x4 

- 

4x3x3 

Swsf 

staggered 

MWwi  Mdw/dz 

- 

- 

3x3x4 

MUw)  Mdu/dz 

1x1x6 

2x2x2 

2x1x1 

S6n 

staggered 

•Mwwy  M. dw/dz 

2x1x1 

1x2x1 

1x1x6 

MUwi  Mdu/dz 

3x3x4 

4x4x4 

4x3x3 

Swn 

staggered 

M^WW  >  Mdw/dz 

4x3x3 

3x4x3 

3x3x4 

39 


TABLE  III:  A-priori  relative  estimation  variation  for  local  estimates,  with  best  global  estimates 


repeated  for  comparison. 


Stencil 

Mfuw 

Mdu/dz 

Mdw/dz 

global  collocated 

0.4430 

0.3530 

0.4660 

0.3137 

global  staggered 

0.4802 

0.4174 

0.4588 

0.3275 

C2 

0.5794 

0.5369 

0.5763 

0.4636 

C4 

0.5185 

0.4620 

0.5160 

0.3744 

C6s 

0.5134 

0.4547 

0.5089 

0.3593 

Cws 

0.5120 

0.4540 

0.5137 

0.3695 

C6n 

0.5134 

0.4518 

0.5089 

0.3593 

Cwn 

0.4590 

0.3638 

0.4998 

0.3289 

S2 

0.5626 

0.5369 

0.5763 

0.4636 

S4 

0.5626 

0.4620 

0.5160 

0.3744 

S6s 

0.5626 

0.4547 

0.5089 

0.3593 

Sws 

0.5099 

0.4575 

0.2638 

0.3695 

S6n 

0.4976 

0.4566 

0.5137 

0.3514 

Swn 

0.2476 

0.4422 

0.4847 

0.3462 
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TABLE  IV:  Coefficients  for  staggered  finite-volume  models,  valid  for  Re\  =  164 


Coefficient 

Standard 

second-order  FV 

Optimal  S2 

Optimal  S4 

A 

v  A"2 

0.314267A-1  +  1.37142i/A~2 

0.491708A-1  +  1.93227i/A~2 

B 

~~ 

- 

-  (0.0941649A-1  -1- 0.26115i/A-2) 

C 

<N 

1 

<1 

0.130001A-1  +  1.24894i/A-2 

0.0963195A-1  +  2.08884z/A-2 

D 

- 

- 

(0.0145042A-1  -  0.19635i/A-2) 

E 

-61/  A-2 

—{2  A  +  AC) 

~(2(A  +  B)  +  4(C  +  D)) 

F 

0.25  A”1 

0.387933  A"1 

0.679643  A"1 

G 

0.5  A"1 

0.270394  A"1 

0.621382  A"1 

H 

- 

- 

0.453439  A"1 

I 

~ 

- 

0.0917196  A"1 

J 

0.25  A"1 

0.268916  A"1 

0.268902  A"1 
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Fig.  1  Stencils  for  1x1x4  and  3x3x4  single  velocity  stencils  are  shown.  The  shaded 
region  identifies  the  face  on  which  flux  is  computed,  (a)  1x1x4  (b)  3x3x4 

Fig.  2  Three-dimensional  energy  spectrum  E(k),  a  -5/3  power  law,  and  the  Nyquist 
wavenumber  associated  with  the  finite- volume  filter  width. 

Fig.  3  Three-dimensional  energy  spectrum  E(k),  filtered  DNS  data  compared  with  opti¬ 
mal  LES  results,  (a)  Models  C2,  C4,  S2,  S4,  S6s  (b)  Models  Cwn,  S6s,  Sws,  S6n,  Swn 
Fig.  4  Third-order  structure  function  Sz(r),  filtered  DNS  data  compared  with  optimal 
LES  results,  (a)  Models  C2,  C4,  S2,  S4,  S6s  (b)  Models  Cwn,  S6s,  Sws,  S6n,  Swn 

Fig.  5  Three-dimensional  energy  spectrum  E(k),  filtered  DNS  data  compared  with  opti¬ 
mal  LES  and  dynamic  Smagorinsky  results. 

Fig.  6  Third-order  structure  function  S3  (r),  filtered  DNS  data  compared  with  optimal 
LES  and  dynamic  Smagorinsky  results. 

Fig.  7  Notation  for  a  one-dimensional  example  problem.  The  state  variables  w  are 
averages  of  some  underlying  conserved  quantity,  whose  flux  across  cell  boundaries  is  denoted 
by  M. 
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FIG.  3:  Zandonade,  Physics  of  Fluids 
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FIG.  6:  Zandonade,  Physics  of  Fluids. 
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FIG.  7:  Zandonade,  Physics  of  Fluids. 
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Abstract 

The  validity  of  the  quasi-normal  approximation,  which  relates  the  fourth-order  velocity 
correlations  to  second-order  velocity  correlations,  is  tested  using  data  obtained  from  direct 
numerical  simulation  of  turbulent  channel  flow  at  ReT  ~  590.  Our  validation  study  indicates 
that  the  quasi-normal  approximation  is  very  accurate  throughout  the  channel  except  for  a 
thin  layer  near  the  wall  {y+  <  50),  where  the  approximation  breaks  down  for  small  separa¬ 
tions  (r+  <  150).  This  study  shows  that  the  quasi-normal  approximation  can  be  used  as  a 
basis  for  development  of  LES  models  using  the  optimal  LES  formulation,  in  wall  bounded 
flows. 

I.  Introduction 

Several  phenomenological  and  statistical  hypotheses  have  been  postulated  to  enhance  our 
understanding  of  turbulent  flows.  One  commonly  used  hypothesis  is  the  Millionshchikov  hy¬ 
pothesis,  which  is  also  known  as  the  quasi-normal  approximation  or  zero-fourth-cummulant 
approximation.  According  to  the  quasi-normal  approximation,  the  fourth-order  velocity 
cummulants  are  zero,  which  allows  us  to  express  the  fourth-order  velocity  correlations  in 
terms  of  second-order  velocity  correlations  (Monin  &  Yaglom1). 

The  quasi-normal  approximation  has  been  used  in  a  variety  of  contexts  in  the  past. 
These  can  be  broadly  categorized  as  (a)  inference  of  turbulence  statistics  from  kinematics 
and  phenomenology  and  (b)  dynamical  evolution  of  turbulence  structure.  In  the  former 
category,  the  quasi-normal  approximation  was  used  by  Batchelor2  to  deduce  the  functional 
behavior  of  pressure  correlations  in  isotropic  turbulence,  Limber3  to  obtain  the  pressure- 
velocity-velocity  correlations  and  Hill  &  Wilczack4  to  comment  on  the  scaling  of  pressure- 
gradient  and  acceleration  in  homogeneous  turbulence.  In  the  latter  category,  the  quasi¬ 
normal  approximation  was  used  by  Proudman  &  Reid5  and  Tatsumi6  to  derive  dynamical 
equations  for  third-order  velocity  correlations  and  Orszag7  in  the  development  of  analytical 
closure  theories  (like  eddy  damped  quasi-normal  markovian  closures). 
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An  experimental  justification  of  the  quasi-normal  approximation  of  fourth-order  veloc¬ 
ity  correlations  was  first  provided  by  Uberoi,8  who  noted  the  conformity  of  his  results  to 
the  quasi-normal  approximation,  except  at  the  small  scales  where  the  approximation  breaks 
down.  However,  the  use  of  quasi-normal  approximation  in  category  (b),  was  questioned 
by  Ogura9  who  noticed  negative  values  in  the  energy  spectrum  after  a  finite  time  during 
the  dynamical  evolution,  a  clear  violation  of  realizability  constraint.  Kraichnan10  argued 
that  the  use  of  quasi-normal  approximation  in  the  Proudman  &  Reid5  approach  was  in¬ 
consistent  with  the  equations  of  motion  and  that  the  inconsistency  results  in  a  violation  of 
energy  conservation.  Orszag7  also  noted  the  pitfalls  of  the  quasi-normal  approximation  in 
dynamical  evolution  equations  and  attributed  the  failures  to  the  dynamical  properties  of  the 
quasi-normal  approximation,  which  include  an  improper  representation  of  relaxation  effects 
resulting  in  a  violation  of  realizability.  Further  attempts  to  correct  these  shortcomings  led 
to  the  development  of  the  eddy-damped  quasi-normal  Markovian  type  closures. 

Most  of  the  objections  to  the  quasi-normal  approximation  are  in  the  context  of  category 
(b)  usages.  However,  inspite  of  the  experimental  justification  of  quasi-normal  approxima¬ 
tion  (Uberoi8),  there  have  been  some  negatives  in  category  (a)  as  well.  It  was  noted  in  Hill 
&  Wilczack4  and  Hill  &  Thoroddsen,11  that  the  quasi-normal  approximation  (although  a 
good  approximation  at  low  Reynolds  numbers)  underestimates  the  pressure  gradient  or  ac¬ 
celeration  correlation  functions  (and  hence  their  variances).  Since  these  correlations  involve 
integral  relations  over  both  small  and  large  scales,  one  should  attribute  the  underestimation 
of  these  correlations  to  the  failure  of  quasi-normal  approximation  at  the  small  scales  where 
intermittency  (non- Gaussian)  effects  are  predominant.  If  one  were  to  consider  quasi-normal 
approximation  for  physical  quantities  that  do  not  depend  on  small  scales  (i.e.  for  separations 
r  »  rj)  one  can  expect  the  approximation  to  be  reasonable. 

In  this  paper,  we  will  examine  the  validity  of  this  proposition  in  a  turbulent  channel  flow 
and  quantify  the  errors  due  to  the  quasi-normal  approximation  using  channel  flow  DNS  data 
at  ReT  ~  590.  The  motivation  for  doing  this  arises  in  a  different  context  than  has  been  cited 
in  the  literature.  Our  objective  is  to  develop  an  optimal  large  eddy  simulation  (LES)  model 
to  simulate  high  Reynolds  number  flow.  Such  a  model  involves  stochastic  estimation  of  the 
convective  acceleration  (or  momentum  flux),  which  requires  the  knowledge  of  fourth-order 
correlations  as  inputs.  These  can  then  be  approximated  in  terms  of  second-order  correlations 
using  the  quasi-normal  hypothesis.  The  channel  flow  was  selected  for  this  study  to  allow 
evaluation  of  the  approximation  of  the  quasi-normal  approximation  in  a  wall  bounded  flow, 
particularly  in  the  log-  or  overlap  region.  The  limits  of  applicability  of  the  approximation 
in  nearly  homogeneous  flows  are  well  understood  as  described  above;  these  limits  need  to 
be  established  for  wall  bounded  flows,  where  models  for  near  wall  turbulence  are  one  of  the 
primary  challenges  in  LES  modeling. 

The  formulation  of  the  optimal  LES  approach,  is  based  on  the  existence  of  an  ideal  LES 
evolution.  The  large  scales  of  turbulence  to  be  simulated  are  defined  through  a  spatial  filter. 
The  choice  of  the  spatial  filter  determines  the  (non-invertible)  mapping  of  the  infinite  dimen¬ 
sional  state  space  of  turbulent  Navier-Stokes  solutions  to  a  finite-dimensional  representation 
of  the  large  scale  field.  For  any  given  spatial  filter,  the  best  deterministic  estimate  (in  the 
mean-square  sense)  of  the  evolution  of  the  large  scale  field,  is  the  ideal  evolution  (Adrian,12 
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Adrian  et  al .,13  Langford  &  Moser14)  which  is  given  by 


dw 

dt 


(1) 


where  u  represents  the  turbulent  field,  w  denotes  the  large  scale  LES  field,  and  •  denotes  the 
filtering  operator.  Such  an  evolution  is  guaranteed  to  produce  the  correct  single-time,  multi¬ 
point  statistics  (Langford  and  Moser,14  Pope15).  Owing  to  the  huge  amount  of  statistical 
information  embedded  in  the  above  conditional  average,  it  is  difficult  to  obtain  it  directly 
and  hence  it  is  formally  approximated  using  the  optimal  LES  approach,  which  is  based  on 
stochastic  estimation.  12,16>13 

Earlier,  this  approach  was  used  for  developing  optimal  sub-grid  models  in  isotropic  turbu¬ 
lence  (Langford  &  Moser14)  and  turbulent  channel  flows  (Volker  et  al.17),  based  on  spectral 
filters.  More  recently,  Langford18  and  Zandonade  et  al.19  developed  a  finite- volume  opti¬ 
mal  LES  technique,  wherein  discretization  and  sub-grid  effects  are  treated  simultaneously. 
Besides  attempting  to  address  the  limitations  (e.g.  filter  definition,  treatment  of  numerical 
discretization  and  near-wall  modeling)  of  other  LES  models,  the  optimal  LES  models  have 
also  led  to  very  accurate  simulations  and  have  the  potential  of  being  applicable  to  the  the 
simulation  of  flows  with  complex  geometries. 

The  optimal  LES  models  require  information  regarding  multi-point  velocity  correlations 
as  inputs.  Although,  these  correlations  can  be  obtained  from  DNS,  it  is  necessary  to  eliminate 
the  need  for  DNS  data,  if  optimal  models  are  to  be  practically  useful.  We  are  thus  led  to 
replace  DNS  statistical  data  with  theory  based  models  of  the  correlations.  For  instance, 
among  the  inputs  needed  for  optimal  LES  models  are  fourth-order  velocity  correlations, 
which  may  be  approximated  using  the  quasi-normal  hypothesis.  The  following  sections  of 
this  paper  explore  the  validity  of  the  quasi-normal  approximation,  within  the  context  of 
optimal  LES  modeling  of  wall-bounded  turbulent  flows. 


II.  Background 

In  the  context  of  our  optimal  LES  models,  we  need  to  determine  an  optimal  quadratic 
(stochastic)  estimate  of  the  convective  acceleration  (or  flux),  which  requires  the  following 
multi-point  correlations  (see  Sec.  2.2.4  in  Zandonade  et  al.19)  as  inputs,  for  any  linear  filter: 

<ui(x)u*(x')) ,  (iti(x)uj(x')ufc(x")),  {ui (x)uj (x')uk (x")ui (x"') )  .  (2) 

Note  that  the  correlations  in  Sec.  2.2.4  in  Zandonade  et  al.19  are  actually  averaged  over  cell 
volumes/faces,  since  an  optimal  convective  flux  term  (averaged  over  a  cell-face)  is  estimated 
in  terms  of  volume  averaged  cell- velocities  (up  to  quadratic  term).  The  volume  averaged 
quantities/correlations  appearing  in  Zandonade  et  al.19  are  a  consequence  of  the  filtering 
operation  inherent  in  the  finite-volume  formulation,  where  the  filter  width  is  related  to  the 
grid  spacing.  Analogous  to  the  finite-volume  optimal  LES  formulation  of  Zandonade  et 
al.,19 one  can  also  attempt  to  construct  optimal  LES  models  based  purely  on  finite  sampling 
of  the  real  velocity  field  (without  the  use  of  a  conventional  filter).  Such  models,  which 
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represent  the  “finite-difference”  analogue  of  finite-volume  optimal  LES,  would  also  require 
correlations  given  in  (2),  as  inputs. 

In  order  to  reduce  the  information  required  as  inputs  (i.e.  in  (2))  for  obtaining  optimal 
quadratic  estimates,  we  use  the  quasi-normal  hypothesis  and  express  the  fourth-order  velocity 
correlations  in  terms  of  second-order  velocity  correlations.  The  most  general  form  of  the 
fourth-order  velocity  correlation  is  given  by 

Qi,j,k,iix » x'>  x"> x"')  =  (x)uj (x')uk (x")ui (x'")) ,  (3) 

which  can  be  approximated  using  the  quasi-normal  hypothesis  as, 

+  Rii(x,x,")Rjk(x',x")> 

(4) 

where  Rij(x,x')  =  (ui(x)uj(xr))  denotes  the  two-point  second-order  correlation.  Note  that 
the  quasi-normal  hypothesis  states  that  the  fourth  order  moments  can  be  determined  as  if  the 
underlying  probability  distribution  were  Gaussian.  However  the  quasi-normal  approximation 
makes  no  statement  on  the  behavior  of  third-order  moments  or  the  underlying  velocity 
distribution. 

To  address  the  modeling  issues  of  wall-bounded  turbulent  flows,  we  use  fully  developed 
turbulent  channel  flow  data  at  ReT  «  590  (Moser  et  al.20)  as  a  test  case  to  quantitatively  val¬ 
idate  the  quasi-normal  hypothesis.  The  data  was  obtained  from  direct  numerical  simulation 
using  periodic  boundary  conditions  in  the  the  streamwise,  (x)  and  spanwise,  (z)  directions, 
while  no  slip  conditions  were  used  on  the  two  parallel  walls.  A  Fourier  spectral  represen¬ 
tation  was  used  in  the  streamwise  and  spanwise  directions  and  a  Chebyshev  representation 
was  used  in  the  wall-normal,  ( y )  direction. 

In  order  to  completely  test  the  quasi-normal  approximation  in  the  most  general  case,  one 
should  compare  the  left  and  right  hand  sides  of  Eq.  4.  However,  this  amounts  to  evaluating  81 
components  of  the  four-point  fourth-order  tensor  at  each  point  in  12-dimensional  space.  This 
evaluation  is  not  practical  for  inhomogeneous  flows,  inspite  of  the  presence  of  symmetries. 

Noting  the  statistical  homogeneities  in  the  streamwise  and  spanwise  directions,  we  restrict 
our  attention  to  a  class  of  two-point  fourth  order  velocity  correlations  to  make  the  study 
tractable.  This  is  a  degenerate  case  of  Eq.  3,  given  by 

QijAx> x0  s  Qi,j,k,i(x> x>  x'>  x0  =  («i(xK(xWx>f  (x')> .  (5) 

Further,  if  the  points  x  and  x'  have  the  same  wall-normal  coordinate  (i.e.  lie  in  a  plane 
parallel  to  the  walls),  we  have 

Qij,ki (x,  x')  =  Qijtki (0,  x'  -  x)  =  QijM  (r) ,  (6) 


Qi,j,k,i(x ,  x'>  x"> x'")  ~  tfii(x>  x,)i?w(x",  x"')  +  Rik{x,  x")Rji(x', x1") 


due  to  statistical  homogeneity  in  the  streamwise  and  spanwise  directions,  where  r  =  x'  —  x. 
Henceforth,  we  will  test  the  quasi-normal  approximation  of  Qij,ki{ r)  for  separation  vectors  r 
parallel  to  the  wall  and  (a)  along  the  streamwise  direction  or  (b)  along  the  spanwise  direction. 
In  order  to  quantify  the  error  due  to  the  quasi-normal  approximation,  we  define  an  error 


measure1,  r),  as, 


fe*(r)  = 


Qij,ki(r)  ~  QNA 


V«,w(y+.r)  may  be  a  better  definition 
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where  QNA  denotes  the  quasi-normal  approximation  of  the  two-point  fourth-order  velocity 
correlation  tensor,  according  to  Eq.  (4)  and  L( r)  =  [<5y,fci(r)  QijMi*)}1^  1S  a  con_ 

traction  (summation  convention  assumed)  or  invariant  of  the  two-point  fourth-order  velocity 
correlation  tensor. 

The  error  measure,  r),  as  defined  in  Eq.  (7),  quantifies  the  error  due  to  the  quasi¬ 
normal  approximation,  associated  with  each  component  of  the  two-point  fourth-order  veloc¬ 
ity  correlation  tensor,  Qtj,w(r)-  Another  error  measure  which  quantifies  the  error  over  all 
the  components  of  Qij,ki( r),  can  be  defined  using  the  same  contraction  of  <£pg,rs(r)  as 

tf(r)  =  \(}>pq,Ts{r)<i>pq,rs{*)}112  •  (8) 

This  error  measure  is  also  relevant  to  the  linear  algebra  of  the  stochastic  estimation  proce¬ 
dure.  The  relative  error  in  the  estimation  coefficients  due  to  the  quasi-normal  approximation 
is  bounded  by  the  condition  number  of  the  linear  system  times  T(r). 

For  separation  vectors,  normal  to  the  wall,  the  fourth-order  velocity  correlations  defined 
in  Eq.  (5)  are  functions  of  two  wall-normal  coordinates,  owing  to  lack  of  homogeneity  in  the 
wall-normal  direction.  In  this  case,  the  fourth-order  velocity  correlations  are  given  by, 

x')  =  Qi,ki(v>  (9) 

where  x  =  ( x,y,z )  and  x'  =  (x,y',z).  Similarly,  as  in  Eq.  (7)  and  Eq.  (8),  we  can  define 
ttjkiM)  and  to  quantify  the  errors  associated  with  the  quasi-normal  approxi¬ 

mation,  for  separation  vectors  perpendicular  to  the  wall. 

In  the  next  section,  the  nature  of  these  errors  will  be  explored  in  the  ReT  w  590  channel 
flow. 

III.  Results 

The  normalized  error  measure  for  wall  normal  velocity  fluctuations,  i.e.  <£22,22 (r),  as  a 
function  of  the  streamwise  separation  (in  wall  units)  and  distance  from  the  wall  ( y+ ,  along  the 
vertical  coordinate)  is  shown  in  Figure  la.  The  separation  vectors,  r,  in  this  case  are  along 
the  streamwise  direction.  The  behavior  of  the  error  as  a  function  of  the  spanwise  coordinate 
is  shown  in  Fig.  lb,  which  is  similar  to  Fig.  la,  except  that  the  separation  vectors  are  chosen 
along  the  spanwise  direction.  The  errors  are  clearly  small  in  these  figures,  indicating  that  the 
quasi-normal  approximation  is  very  good  for  wall  normal  fluctuations  with  separations  along 
the  spanwise  coordinate.  The  normalized  error  for  spanwise  fluctuations  (<£33,33 (*))  is  also 
very  small  as  shown  in  figure  2.  Although,  the  data  of  Moser  et  al.20  indicate  large  values 
of  flatness  factors  of  wall-normal  (and  spanwise)  velocity  fluctuations  near  the  wall,  which 
significantly  differ  from  the  quasi-normal  value,  the  error  measure  that  is  relevant  for  our 
validation  of  the  quasi-normal  approximation  for  wall-normal  (and  spanwise)  fluctuations 
within  the  context  of  optimal  LES  is  <£22,22  (r). 

The  error  measure  <£n,ii(r)  for  streamwise  and  spanwise  separations  (similar  to  Figs,  la 
and  lb)  are  shown  in  figures  3a  and  3b  respectively.  The  contour  levels  in  Fig.  3a  indicate 
that  the  error  is  small  everywhere  except  for  regions  near  the  wall  ( y+  <  30)  and  for  small 
streamwise  separations  (Aa;+  <  150),  where  the  error  is  about  20%  .  A  similar  observation 
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can  also  be  made  from  Fig.  3b,  where  the  normalized  error  is  small  everywhere  except  in 
the  near- wall  region  (y+  <  50)  and  for  small  spanwise  separations  (Az+  <  100).  The 
deviation  from  quasi-normality  in  a  region  very  near  the  wall  should  be  expected  from 
the  intermittent  nature  of  the  near  wall  flow.  Although  the  quasi-normal  approximation 
breaks  down  for  small  streamwise  and  spanwise  separations,  our  model  inputs  require  fourth 
correlations  approximated  by  quasi-normal  approximation  only  for  streamwise  separations 
that  are  greater  than  those  in  the  breakdown  region.  Typical  channel  flow  optimal  LES 
simulations  (Volker  et  al.  2002)  would  require  correlations  with  separations  of  order  Ax+  > 
300  and  A z+  >  150. 

To  confirm  the  good  accuracy  of  the  quasi-normal  approximation  described  above,  a 
number  of  other  components  of  the  fourth  order  correlation  tensor  were  checked  and  found 
to  have  small  errors  due  to  quasi-normality.  Further,  the  overall  error  in  the  representation 
of  the  fourth  order  correlation  tensor  as  measured  by  d'(r),  is  shown  in  figure  4.  Note  that 
the  contours  of  T  shown  in  figure  4  are  similar  to  those  for  <j)u,u  in  figure  3.  Thus  the 
error  in  the  quasi-normal  approximation  in  this  flow  is  dominated  by  the  streamwise  velocity 
fluctuation  (0n,n). 

For  separation  vectors  in  the  wall-normal  direction,  the  errors  022, 22(2/)  v')i  033,33 (?/>  v') » 
Wi  niv^y')  and  '&±(y,y')  are  Shown  in  figure  5.  With  the  exception  of  the  region  very  near 
the  wall  (y+  <  40,  yl+  <  40  ),  the  errors  are  again  very  small  and  the  overall  error  as 
measured  by  '&±(y,y')  is  dominated  by  the  streamwise  fluctuations  (0n,ii). 

IV.  Discussion 

In  this  paper,  we  have  used  data  obtained  from  direct  numerical  simulation  of  turbulent 
channel  flow  at  ReT  «  590,  to  test  the  validity  of  the  quasi-normal  approximation  of  fourth- 
order  velocity  correlations  (see  Eq.  (4)).  Unlike  previous  studies  (mentioned  in  Sec.  1),  our 
motivation  arises  in  the  context  of  stochastic  estimation  of  the  convective  acceleration  (or 
momentum  flux)  term,  in  order  to  develop  optimal  LES  models. 

The  error  due  to  the  quasi-normal  approximation,  of  a  class  of  two-point  fourth  order 
correlations,  was  quantified  using  a  component- wise  error  measure  r)  and  an  invari¬ 
ant  error  measure  T(r)  for  separation  vectors  in  the  streamwise,  spanwise  and  wall-normal 
directions.  Our  results  indicate  that  the  quasi-normal  approximation  provides  an  excellent 
representation  of  the  fourth-order  correlation  tensor,  except  in  the  region  very  near  the 
wall  (y+  <  50)  where  there  are  significant  errors  in  the  approximation  for  small  separation 
(r+  <  150).  In  this  region,  the  error  in  the  representation  is  dominated  by  streamwise 
velocity  fluctuations. 

Note  that  the  error  measures  used  here  are  appropriate  for  evaluating  the  representation 
of  the  correlation  tensor  as  a  whole.  Because  the  errors  are  normalized  by  a  norm  of  the 
tensor,  small  values  of  a  component  error  measure  (e.g.  022, 22)  does  not  imply  that  the 
related  error  in  that  individual  component  is  necessarily  small.  For  example,  022,22  is  small 
near  the  wall,  but  it  is  known  that  the  flatness  factor  of  the  wall-normal  fluctuation  u2) 
denoted  by  Pi{u2)  =  (u\)/{u2)2,  is  far  from  the  Gaussian  value  of  3  (according  to  the  data  of 
Moser  et  al.20).  Thus  {u\)  is  by  itself  not  well  represented  by  the  quasi-normal  approximation 
of  3(1*2)  at  the  wall.  The  error  measure,  022,22  is  small  in  this  case  because  near  the  wall 
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fluctuations  of  u2  are  insignificant  compared  to  those  of  u\  (i.e.  streamwise  fluctuation), 
since  U2  goes  to  zero  like  y2  near  the  wall.  Thus  near  the  wall,  U\  fluctuations  dominate  the 
norm  of  the  fourth  order  correlation  tensor.  These  error  measures  are  relevant  to  our  use 
of  the  quasi-normal  approximation  to  model  statistical  quantities  for  input  to  optimal  LES. 
The  stochastic  estimates  performed  in  optimal  LES  use  fourth  order  correlation  components 
(or  their  integrals)  as  elements  of  the  right  hand  side  vector  in  a  linear  algebraic  equation 
for  the  estimation  kernels.  In  this  context,  the  L 2  norm  of  the  error  in  the  quasi- normal 
approximation  will  be  directly  related  to  the  L2  norm  in  the  resulting  estimation  kernel, 
which  in  turn  is  directly  related  to  the  error  due  to  this  approximation  in  the  LES  model. 

The  validity  of  quasi-normal  approximation  in  wall  bounded  turbulence,  except  for  a 
very  thin  layer  near  the  wall,  provides  one  of  several  theoretical  modeling  tools  that  will 
be  required  to  make  optimal  LES  a  viable,  practical  modeling  paradigm  for  wall  bounded 
turbulence. 
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Figure  4a:  Contours  of  the  invaria: 
distance  from  the  wall,  as  in  figure 
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Figure  5d:  Contours  of  the  invariant, 
direction 


